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ABSTRACT 

This  report  documents  a  program  for  the  numerical  solution  of 
large  sparse  systems  of  algebraic  and  implicitly  defined  stiff  differen¬ 
tial  equations.  The  principal  use  is  intended  to  be  the  solution  of 
differential  equations  arising  from  time  dependent  partial  differential 
equations  when  the  finite  element  method  is  used  to  discretize  the 
space  domain.  The  use  of  compact  matrix  storage  techniques  and  iteration 
for  the  solution  of  the  quasi-Newton  iterates  in  Gear’s  method  makes 
the  program  extremely  efficient  both  in  terms  of  storage  requirements 
and  execution  times. 
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1,0  Introduction 


This  report  documents  a  program  for  the  solution  of  algebraic  and 
implicitly  defined  stiff  differential  equations.  We  were  particularly 
interested  in  solving  very  large  systems  of  differential  equations 
arising  from  partial  differential  equations  where  the  finite  element 
method  has  been  applied  in  the  space  variables. 

Our  original  goal  was  to  use  a  compact  storage  scheme  for  the 
large  matrices  involved  and  to  use  iteration  to  solve  the  linear 
algebraic  systems  which  occur.  However,  the  resulting  program  is  easily 
adapted  to  different  applications  through  user  modifications  accomplished 
by  replacement  of  one  or  two  relatively  simple  subroutines.  Thus  the 
program  is  a  powerful  one  which  can  be  used  in  a  variety  of  applications. 
Four  examples,  illustrating  different  matrix  storage  techniques  and 
different  linear  equation  solvers  are  given  in  the  appendix.  Other 
storage  schemes  and  solution  methods,  e.g.  symmetric  Jacobian  with 
Cholesky  decomposition,  are  relatively  simple  to  implement. 

In  Section  2  a  brief  review  of  the  integration  scheme  is  given. 

In  Section  3  a  discussion  of  differences  between  this  program  and  the 
one  from  which  it  was  adapted  is  given.  Section  4  is  devoted  to  a 
discussion  of  information  concerning  the  use  of  this  package. 

2.0  Theoretical  Background 

Consider  the  system  of  N  =  m  +  il  differential  and  algebraic 
equations  in  •  •  •  » 

(1)  F(y,y,t)  +  P(t)  V  =  0  , 


1 


with  all  or  some  of  the  initial  values  y^^Ctg)  , . .  .  ,y^(tQ)  , . . .  , 

v^(tg)  specified.  Enough  of  the  above  values  must  be  given  in  order 

to  determine  the  remaining  values  and  initial  values  for  any  of  the 

•  • 

derivatives,  appears  in  equation  (1).  In  equation  (1) , 

P(t)  is  an  N  X  matrix,  F  is  a  vector  with  N  components,  and 
V  is  a  vector. 

The  program  documented  here  is  a  modification  of  a  program  due  to 
Brown  and  Gear  [1].  The  method  of  solution  is  a  modification  of  Gear’s 
method  for  stiff  differential  equations  [2].  The  application  to  differ¬ 
ential  algebraic  systems  was  given  by  Gear  [3].  We  will  briefly  describe 
the  method  here  for  completeness,  and  refer  the  reader  to  the  references 
for  more  details. 

Suppose  that  approximate  solution  values  are  known  at  a  number  of 
equally  spaced  points,  ^n-l’^n  2*****^n-k*  represented  by 

y^’^  ^\...,y^^  respectively.  Let  ^(^^-1^  represented  by 

v(n  1)  ^  ^  backward  differentiation  formula  gives 

,  •  (n)  1  f  (n)  (n-1)  (n-k)  v 

h  y^  =  (ttg  y  +  y  +  .  . .  +  otj^  y  )  y 

where  h  =  coefficients  and  6q  are  from  Gear 

[2] ,  p.  217.  Substitution  of  this  into  (1)  gives 


(2)  ,  t^)  +  P(t^)  =  0 


6oh 


n 


as  the  equation  which  y^'^^  and  must  satisfy.  In  equation  (2)  , 


V  1  /  (n-1)  .  .  (n-k). 

+  ...  +  aky  )  • 
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In  general,  equation  (2)  represents  a  system  of  nonlinear  equations  for 
and  •  The  method  used  for  solving  this  system  of  algebraic 

equations  is  a  variant  of  Newton’s  method.  The  initial  guess,  , 

is  obtained  by  polynomial  extrapolation  using  a  Hermite  polynomial 
interpolating  the  known  values  y^’^  ...,  y^’^  .  Thus 

the  predicted  values  has  the  form 


(3) 


=  h 


6,  y 


(n-l) 


+  a,  y 


(n-l) 


a  y 
n-k  ■' 


(n-k) 


The  application  of  Newton's  Method  to  equation  (2)  then  yields  the 
corrector  equation 


(4)  J  • 


(n)  ,  i+1  _  (n)  ,  i 

=  -  F(y^"^’^,  t^)+P(t^)  , 

^(n) ,i+l  _  ^(n) ,i  /  ^0 


where  J  is  the  Jacobian  matrix. 


(5) 


Gear  shows  that  the  initial  guess  for  is  not  important, 

and  thus  is  used.  Up  to  three  iterations  are  performed  on  the 

corrector  equation.  The  matrix  J  is  not  evaluated  at  each  iteration, 
nor  even  at  each  timestep.  J  is  evaluated  whenever  (i)  the  timestep 
or  order  is  changed,  or  (ii)  the  corrector  iteration  fails  to  converge 
in  three  iterations. 

If  the  corrector  iteration  fails  to  converge,  the  J  matrix  is 
evaluated,  unless  it  had  already  been  evaluated  at  the  current  time. 
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If  it  has  been  evaluated  at  the  current  time,  the  timestep  is  reduced 
by  a  factor  of  4,  In  either  case  the  step  is  then  retired. 

If  the  corrector  iteration  converges,  the  local  error  is  estimated, 
based  on  the  fact  that  the  local  error  is  approximately  proportional  to 
the  difference  between  the  predicted  and  corrected  values  of  y^^^  • 

For  this  purpose  a  relative  error  tolerance  is  used  for  large  solutions 
and  an  absolute  error  tolerance  for  small  solutions.  The  root-mean- 
square  norm  (Euclidean  norm  divided  by  the  square  root  of  the  number  of 
components)  is  used  for  the  vector  with  components  e^/ymax^  ,  where 
e^  is  the  estimated  local  error  in  ymax^  =  max  (|yf^^|,  1)  . 

If  the  error  is  larger  than  that  specified  by  the  user,  an  acceptable 
timestep  is  estimated  for  the  current  order  or  order  one  lower,  and 
the  step  repeated.  Up  to  three  such  failures  are  permitted,  after  which 
an  attempt  is  made  to  start  over  with  a  first  order  method. 

When  using  a  method  of  order  q  ,  the  program  takes  at  least 
q  +  1  steps  before  changing  the  timestep.  Changes  in  timestep  are 
preceeded  by  calculation  of  the  predicted  timestep  at  current  order  and 
order  one  higher  and  one  lower.  If  the  timestep  can  be  increased  by  more 
than  10%,  the  order  corresponding  to  the  largest  possible  timestep  is 
used.  If  the  timestep  cannot  be  increased  by  at  least  10%,  the  current 
order  and  stepsize  are  retained  for  at  least  10  more  steps. 

After  each  step  a  test  is  made  to  determine  whether  time  has 
advanced  to  or  beyond  the  input  end  time.  Control  is  returned  to  the 
calling  program  if  it  has. 
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At  the  initial  call,  no  history  of  the  solution  is  available,  so 
the  program  must  begin  with  a  first  order  method,  taking  two  such  steps 
in  accordance  with  the  above  description.  The  times tep  must  be  suitably 
small,  again  in  accordance  with  the  above.  At  the  point  the  program  can 
begin  to  increase  the  order  of  the  method  and  the  times tep.  Because 
the  Jacobian  matrix  must  be  generated  whenever  the  timestep  is  changed, 
it  is  not  efficient  to  try  too  large  a  timestep  initially.  Because  the 
program  rapidly  finds  the  best  order  and  timestep,  it  is  relatively  cheap 
to  underestimate  the  initial  timestep  compared  to  the  cost  of  overestimating 
it. 


3,0  Differences  compared  with  DFASUB 

The  principal  differences  between  the  SDESOL/LDASUB  package  and 
DFASUB,  and  the  reasons  for  incorporating  them  are  as  follows. 

(i)  The  main  goal  of  this  revision  was  to  generate  a  program  which 
could  solve  very  large  sparse  systems  of  differential  equations  efficiently, 
both  in  terms  of  storage  requirements  and  execution  time.  We  are 
particularly  interested  in  the  solution  of  ordinary  differential  equa¬ 
tions  arising  from  partial  differential  equations  where  the  finite 
element  method  has  been  used  to  discretize  the  problem  in  space. 

Large  sparse  problems  require  at  least  a  different  system  of  stor¬ 
ing  the  Jacobian  matrix  and  possibly  the  use  of  iteration  to  solve  for 
the  quasi-Newton  iterates  in  equation  (2.4).  Two  such  subroutine 
packages,  to  be  used  with  the  basic  subroutines,  have  been  provided. 

Another  package  using  standard  elimination  techniques  is  also  provided 
and  is  convenient  for  smaller  systems  of  equations.  Use  of  any  of  these 
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options  requires  the  user  to  supply  a  subroutine  to  evaluate  his  form 
of  the  equations  (1),  and  for  efficiency,  a  subroutine  to  explicitly 
evaluate  the  Jacobian •  A  subroutine  to  approximate  a  full  Jacobian 
through  numerical  differencing  has  been  provided.  With  the  exception 
of  a  minor  correction,  this  is  the  same  as  given  in  [1].  While  use  of 
this  routine  is  convenient,  it  is  inefficient  and  should  be  avoided 
for  large  systems.  It  is  anticipated  that  the  user  can  provide  his  own 
subroutine  package,  using  his  own  storage  scheme  for  the  Jacobian,  and 
with  a  suitable  equation  solver  for  the  Newton  iterates.  There  should 
be  no  need  to  disturb  the  basic  package  which  carries  out  the  time 
integration. 

(ii)  For  user  convenience,  without  a  major  rewrite  of  DFASUB,  a  driver 
routine,  SDESOL,  to  be  referenced  by  the  user  and  which  then  communicates 
with  LDASUB  was  written.  The  chief  function  of  SDESOL  is  to  set  up  a 
number  of  references  to  work  storage  areas  required  by  LDASUB.  In 
addition,  some  testing  of  parameters  is  accomplished,  and  a  subroutine 

to  calculate  initial  values  of  derivatives  is  called. 

(iii)  A  subroutine  to  calculate  initial  values  of  derivatives,  DERVAL, 

has  been  provided.  The  routine  provided  requires  that  the  first  m 
9F 

rows  of  —  be  nonsingular,  which  does  not  need  to  be  true  in  the 

9y 

general  case.  For  this  reason,  and  others  discussed  in  Section  4,  the 
user  may  need  to  provide  either  his  own  version  of  DERVAL  to  evaluate 
the  derivatives  initially,  or  else  he  may  supply  initial  values  and  a 
dummy  version  of  DERVAL. 

(iv)  Other  changes  made  in  generating  LDASUB  from  DFASUB  were  to  simplify 
the  code  for  the  particular  type  of  problem  we  wish  to  solve,  while 
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others  were  to  enhance  the  usability  of  the  code.  Some  errors  were 
also  corrected,  notably  two  errors  in  coefficients  for  the  fifth  and 
sixth  order  methods.  DFASUB  had  the  capability  of  computing  various 
elements  of  the  Jacobian  at  different  times  if  they  had  different 
dependencies,  with  the  possibility  of  inverting  that  part  of  the  matrix 
at  that  time,  if  it  could  be  done.  This  could  result  in  increased 
efficiency  in  certain  problems,  at  some  expense  in  convenience,  but  for 
our  purposes  it  was  not  considered  useful,  and  was  removed.  Therefore 
only  one  call  is  made  to  evaluate  the  Jacobian.  The  Jacobian  was 
evaluated  at  the  beginning  of  each  times tep  in  DFASUB,  but  this  has 
been  eliminated  in  LDASUB,  in  accordance  with  the  description  in  Section 
2.  A  subroutine,  S2,  was  called  from  DFASUB  to  evaluate  time  dependent 
terms  whenever  time  was  changed.  This  is  reasonable,  since  the  function 
evaluation  routine  may  be  called  several  times  at  a  given  value  of  the 
independent  variable.  We  have  removed  this,  preferring  to  test  for  a 
new  time  in  the  routines  where  time  dependent  terms  appear,  then 
evaluating  and  storing  them  internally  to  that  routine  when  necessary. 
This  helps  make  the  function  evaluation  more  self  contained,  as  well. 

In  DFASUB  extra  parameters  in  the  calling  sequence  allowed  the 
user  to  communicate  constants  to  the  function  evaluation  and  Jacobian 
subroutines  via  DFASUB.  We  believe  this  is  inefficient  and  confusing, 
and  removed  this  capability,  preferring  to  communicate  from  the  main 
program  to  these  subroutines  via  Common,  or  possibly  through  multiple 
entry  points. 

The  norm  used  for  error  tests  in  DFASUB  is  the  Euclidean  norm. 

This  has  the  undesirable  property  that  for  large  systems  the  allowable 
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error  criterion  may  be  large.  We  therefore  changed  to  the  root^mean- 

square  norm  in  LDASUB,  which  is  simply  the  Euclidean  norm  divided  by 

the  square  root  of  the  number  of  components.  One  other  change  was 

made  in  the  error  tests.  As  noted  in  Section  2,  the  error  vector  has 

components  e^/ymax^  ,  where  e^  is  the  estimated  local  error  in  the 

i—  variable  y.  ,  and  ymax,  =  max  (|y^^^|,  1)  .  In  DFASUB  the 
^  0_<k_<n  ^ 

maximum  was  taken  only  up  to  the  previous  timestep,  n-1  .  This  change 

was  incorporated  because  some  of  the  problems  in  which  we  were  interested 

began  with  many  components  at  zero,  but  which  very  rapidly  became  large, 

12 

around  10  or  more.  Without  updating  the  value  of  ymax^  ,  the 
size  of  the  timestep  was  artificially  kept  extremely  small  in  order  to 
satisfy  an  unreasonable  error  tolerance.  For  this  reason,  the  maximum 
value  of  the  component  was  updated  before  the  norm  of  the  relative  errors 
was  computed.  For  problems  where  values  range  near  to  one,  the  modifica¬ 
tion  will  result  in  no  appreciable  change  in  performance. 

The  printout  of  counters,  timestep,  time,  and  values  of  the  dependent 
variables  was  made  an  option  through  a  parameter  in  the  calling  sequence. 

An  additional  value  printed  is  the  order  of  the  method  being  used. 

DFASUB  incorporated  the  capability  of  terminating  if  a  certain 
number  of  floating  point  overflows  had  occurred.  This  capability  was 
removed  from  LDASUB. 

The  final  modification  to  the  program  was  the  incorporation  of  a 
restart  capability  without  having  to  begin  again  with  a  first  order 
method.  This  was  accomplished  by  adding  two  entry  points  to  LDASUB. 

One,  LDASAV,  saves  values  internal  to  LDASUB,  returning  them  to  the 

main  program,-  where  they  can  be  saved  for  the  time  at  which  the  calculation 
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is  to  be  resumed.  At  that  time,  another  entry  point,  LDARST,  restores 
those  values  internal  to  LDASUB,  while  other  necessary  values  are 
restored  through  the  calling  sequence. 

4.0  Subroutine  Descriptions 

The  description  of  subroutines  is  divided  into  two  subsections. 

The  first  deals  with  the  basic  integration  routine  and  other  subroutines 
which  make  up  the  core  package,  and  which  should  not  be  modified  by  the 
casual  user.  The  second  deals  with  a  set  of  supporting  subroutines,  at 
least  one  of  which  must  be  supplied  by  the  user  since  it  defines  his 
system  of  equations.  The  others  may  be  usable  in  the  form  given  in 
one  of  the  examples,  or  can  be  defined  by  the  user  to  accomplish  his 
desired  implementation. 

4. 1  Basic  Subroutine  Package 

4.1.1  Subroutine  SDESQL 

This  routine  is  the  only  one  which  needs  to  be  referenced  by  most 
users.  It  serves  as  a  driver  for  the  integration  routine,  LDASUB. 

SDESOL  has  a  simpler  calling  sequence  than  LDASUB  and  relieves  the 
user  of  having  to  set  up  a  number  of  auxiliary  storage  arrays.  In 
addition,  the  routine  calls  DERVAL  to  calculate  the  values  of  the 
derivatives  on  the  initial  call. 

The  calling  sequence  is 

CALL  SDESOL (Y,YL,T, TEND, NY, NL,M,JSKF,MAXDER,IPRT,H,HMIN,HMAX,RMSEPS,W) 
where  the  parameters  are  defined  as  follows. 


9 


Y  -  Input  and  output.  An  array  dimensioned  (7, NY).  On  the  initial 
call  this  array  contains  the  initial  values  of  the  dependent 


stored  in  Y(j+l,i)  ,  Here  h  is  the  current  stepsize. 

These  values  must  not  be  changed  between  returns  to  the 
calling  program  and  subsequent  entries  to  SDESOL.  It  is 
possible  to  interpolate  for  values  of  the  dependent  variables 
at  times  other  than  those  calculated  by  using  the  formula 


,  where  q  is  the  order  formula 


j=0 


currently  in  use,  and  is  obtained  as  q  =  |JSKF/10|  . 

YL  -  input  and  output.  Array  of  linear  variables,  v^,  i=l>  •••,  ^ 
The  user  supplies  initial  values  for  these  variables,  and 
during  execution  it  contains  current  values  of  the  linear 
variables. 

T  -  input  and  output.  The  user  supplies  the  initial  time,  which 
is  updated  to  current  time  during  execution. 

TEND  -input*  Time  at  which  the  integration  is  to  end.  This  is  the 
only  parameter  normally  changed  by  the  user  between  succes¬ 
sive  entries  to  SDESOL. 

NY  -  input.  The  number  of  differential  and  nonlinear  variables,  m 

NL  -  input.  The  number  of  linear  variables,  H  .  This  may  be  zero. 

M  -  input.  The  number  of  variables  to  be  included  in  the  local 

error  test.  The  error  test  will  be  performed  for  the  variables 
y^,  i=l,  M  .  The  M  used  is  no  greater  than  NY  . 
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JSKF  -input  and  output.  An  indicator:  on  input, 

JSKF  =  0  indicates  that  this  is  the  initial  call  to  SDESOL. 
Initial  values  of  the  derivatives  are  calculated  and  auxiliary 
storage  references  are  set  up.  This  also  indicates  to 
subroutine  LDASUB  to  initialize  parameters  and  begin  with 
a  first  order  integration  method. 

JSKF  >  0  indicates  a  continuation  of  a  previous  call  to  SDESOL 
JSKF  =  -  1  indicates  a  restart  call.  This  is  discussed 
further  in  Section  4.1.2. 

JSKF  <  -  1  may  result  from  the  user  neglecting  to  test  for 
error  returns  from  SDESOL.  Because  of  this  possibility, 
the  run  is  terminated  with  an  appropriate  comment  when 
JSKF  <  “  1  is  input. 

On  output,  JSKF  normally  is  a  two  digit  number,  ±  qp  .  q 
is  the  order  of  the  formula  currently  being  used  for  the 
integration,  p  is  an  indicator  determining  the  type  of 
return.  JSKF  >  0  ,  p  =  1  is  the  normal  return.  Note  that 
SDESOL  may  be  re-entered  to  continue  the  solution  without 
changing  JSKF.  JSKF  <0  is  an  error  return,  with  the  value 
of  p  indicating  the  error,  as  follows, 
p  =  1  error  test  failure  for  H  ^  HMIN 
p  =  3  corrector  failed  to  converge  for  H  >  HMIN 

p  =  4  corrector  failed  to  converge  for  a  first  order  method 

p  =  5  error  return  from  subroutine  NUITSL 

p  =  6  error  return  from  subroutine  DERVAL 
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MAXDER- input .  Maximum  order  method  which  should  be  used.  The 
highest  order  possible  is  six. 

IPRT-  input.  Print  control  indicator. 

0  ,  no  print  from  LDASUB 

>  0  ,  at  each  step,  print  number  of  steps,  number  of  Jacobian 
evaluations , current  order  being  used,  stepsize  for  next  step, 
current  time,  and  current  values  of  the  dependent  variables. 

H  -  input  and  output.  On  initial  call  it  is  an  estimate  of  the 

timestep.  During  execution  it  is  updated  to  the  current  value, 
and  on  return  contains  the  stepsize  to  be  tried  for  the  next 
step.  The  input  value  need  not  be  accurate.  It  is  better  to 
underestimate  than  to  overestimate  the  initial  value.  The 
stepsize  and  order  are  varied  to  meet  the  local  error  tolerance 
specified.  The  user  does  not  normally  change  the  stepsize 
between  entries  to  SDESOL. 

HMIN-  input.  The  minimum  stepsize  to  be  allowed. 

HMAX-  input.  The  maximum  stepsize  to  be  allowed. 

RMSEPS-input .  The  error  test  constant.  The  values  of  the  relative 
local  errors  must  have  root -mean-square  norm  less  than  RMSEPS. 

W  -  an  array  of  auxiliary  storage  required  by  LDASUB.  This  array 
must  contain  a  total  number  of  locations  equal  to  the  sum  of 
(i)  13*NY  +  5*NL  for  arrays  used  in  LDASUB,  (ii)  storage 

for. the  Jacobian  matrix,  and  (iii)  any  locations  used  in 
processing  the  Jacobian,  e.g.,  scratch  storage  used  by  an 
equation  solver. 
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4. 1,2  Subroutine  LDASUB 


This  subroutine  is  the  basic  integration  routine  and  performs  the 
process  in  essentially  the  same  manner  as  subroutine  DFASUB.  A  brief 
description  is  given  in  Section  2  and  differences  between  this  routine 
and  DFASUB  are  outlined  in  Section  3.  Parameters  in  this  routine  in 
which  the  user  may  be  interested  are  stored  in  the  W  array,  an  argument 
of  subroutine  SDESOL, 

YMAX  -  array  of  maximum  magnitudes  of  the  independent  variables, 
y^,  up  to  the  current  time  (or  one,  if  less  than  one). 

This  is  stored  beginning  at  location  7*NY  +  NL  +  1  of 
the  W  array. 

ER  “  This  is  the  array  of  differences  between  the  predicted 
and  corrected  values  of  the  variables,  y^  ,  and  is 
proportional  to  the  estimated  error.  This  array  is 
stored  beginning  in  location  8*NY  +  NL  +  1  of  the  W 
array. 

This  subroutine  incorporates  a  restart  capability.  In  order  to 
restart  from  a  previous  point  without  beginning  again  with  a  first  order 
method,  it  is  necessary  to  have  saved  a  number  of  variables  internal  to 
LDASUB,  and  then  restore  them  before  calling  SDESOL  again.  To  save  the 
internal  parameters,  the  user  calls  subroutine  LDASAV(SAV).  Here  SAV 
is  an  array  of  length  29  in  which  the  values  to  be  saved  will  be  stored. 
In  addition  to  SAV,  the  user  must  also  save  a  number  of  arrays  in  the 
calling  sequence  of  LDASUB,  and  this  is  most  easily  accomplished  by 
saving  the  W  array  in  the  calling  sequence  for  SDESOL.  Once  these 
arrays  have  been  saved,  along  with  the  other  simple  parameters  in  the 
calling  sequence  (  Y  and  YL  need  not  be  saved) ,  the  user  is  free  to 
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use  the  package  to  solve  a  different  problem,  or  to  terminate  the  computer 
run,  to  be  restarted  later. 

At  the  time  the  problem  is  to  be  restarted,  the  user  calls  sub¬ 
routine  LDARST(SAV),  where  SAV  is  the  array  of  values  obtained  previously 
by  calling  LDASAV.  This  restores  internal  values  in  LDASUB.  The  user 
then  calls  SDESOL  with  the  same  simple  parameters  and  the  W  array  as 
before,  except  that  JSKF  =  -  1  and  a  new  end  time,  TEND,  is  provided. 
Restoration  of  values  (including  Y  and  YL  )  in  LDASUB  is  completed 
and  solution  of  the  problem  resumes. 

If  the  user  desires  to  change  the  error  tolerance,  number  of 
variables  in  the  error  test,  or  maximum  order  to  be  used,  the  user 
must  make  a  new  initial  call  to  SDESOL,  that  is,  set  JSKF  =  0  . 

4.1.3  Subroutine  CQPYZ 

This  subroutine  simply  transfers  the  contents  of  one  array  into 
another  array. 

4.2  Supporting  Subroutines 

This  group  of  subroutinesmust ,  at  least  in  part,  be  supplied  by 
the  user.  The  user  must  supply  at  least  one  subroutine,  DIFFUN.  For 
better  efficiency,  the  user  should  supply  a  subroutine,  JACMAT,  to 
explicitly  evaluate  the  Jacobian,  although  a  version  which  approximates 
the  Jacobian  by  numerical  differencing  is  given  in  the  appendix.  To 
take  advantage  of  sparsity  or  other  features  of  his  problem,  the  user 
will  need  to  supply  the  subroutine  NUITSL  to  solve  the  systems  of 
equations  (2.4).  For  certain  problems  the  user  may  have  to  supply 
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subroutine  DERVAL  to  calculate  the  initial  values  of  the  derivatives. 
We  discuss  the  requirements  of  these  subroutines  in  turn. 


4.2.1  Subroutine  DIFFUN 

This  subroutine  simply  evaluates  the  equations  (2.1)  at  a  given 
time  and  yalues  of  y  ,  y  ,  and  V  .  Other  parameters  in  the  function 
definition  must  be  transmitted  from  the  calling  program  via  COMMON  or 
some  other  device,  determined  by  the  user. 

The  calling  sequence  is 

CALL  DIFFUN  (Y,  YL,  T,  HINV,  DY) ,  where  the  parameters  are  defined  as 
follows . 


Y 


YL 


T 

HINV 

DY 


input.  Same  as  in  SDESOL.  This  array  contains  the  current 
values  of  the  variables  y^  and  their  (scaled)  derivatives, 
input.  Same  as  in  SDESOL.  This  array  contains  the  current 
values  of  the  linear  variables, 
input.  Current  time. 

input.  1/h  ,  where  h  is  the  current  stepsize. 
output.  Array  of  function  values. 


4.2.2  Subroutine  JACMAT 

This  subroutine  evaluates  the  Jacobian  matrix  J  ,  equation  (2.5) 
at  the  given  time  and  current  values  of  the  dependent  variables,  order, 
and  stepsize.  A  version  of  JACMAT  which  approximates  J  by  numerical 
differencing  is  given  in  the  appendix.  For  maximum  efficiency,  the 
user  should  supply  the  explicit  representation  of  the  Jacobian.  Because 
the  Jacobian  is  used  to  solve  for  the  quasi-Newton  iterates,  it  is  not 
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necessary  for  the  Jacobian  to  be  exact.  Thus  the  user  should  consider 


the  possibility  of  approximations  which  reduce  the  total  number  of 
computations  in  this  step,  with  due  regard  for  the  fact  that  a  smaller 
timestep  may  be  required  to  obtain  convergence  of  the  corrector  within 
three  iterations. 

The  calling  sequence  for  this  subroutine  is 
CALL  JACMAT  (Y,  YL,  T,  HINV,  A2,  N,  NY,  EPS,  DY,  FI,  PW) ,  where  the 
parameters  are  defined  as  follows. 


Y 

input.  Same  as  in  SDESOL,  Y  contains  the  current  values 

YL 

of  the  variables  y^  and  their  (scaled)  derivatives. 

input.  Same  as  in  SDESOL.  This  array  contains  the  current 

values  of  the  linear  variables. 

T 

input.  Current  time. 

HINV  - 

input.  1/h  ,  where  h  is  the  current  stepsize. 

A2 

input.  The  constant  from  LDASUB. 

N 

input.  Total  number  of  variables. 

NY 

input.  Number  of  differential  and  nonlinear  variables. 

EPS 

input.  Error  constant  from  LDASUB,  ’RMSEPS. 

DY 

input.  Array  of  current  function  values. 

FI 

scratch  array  of  N  locations  available  for  use  by  this 

routine. 

PW 

output.  The  Jacobian  matrix  J  ,  or  an  approximation, 

calculated  in  JACMAT  and  returned  to  calling  program. 

This  matrix  is  used  in  subroutine  NUITSL  and  the  storage 

mode  must  agree  between  the  two  subroutines. 
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4.2.3  Subroutine  NUITSL 

This  subroutine  solves  the  equations  (2.4)  for  the  quasi-Newton 
iterates.  This  subroutine  will  normally  be  supplied  by  the  user, 
although  versions  which  solve  the  system  by  elimination  methods  and 
iterative  methods,  respectively,  are  given  in  the  examples  in  the 
appendix.  This  subroutine  will  often  be  modified  or  replaced  by  the 
user  to  take  advantage  of  sparsity  or  other  features  of  his  problem  in 
connection  with  JACMAT,  of  course. 

The  calling  sequence  for  this  subroutine  is 
CALL  NUITSL  (PW,  DY,  FI,  N,  NY.  EPS,  YMAX,  NEWPW,  KRET) ,  where  the 
parameters  are  defined  as  follows. 

PW  -  input.  The  Jacobian  matrix  J  computed  in  JACMAT. 

DY  “  input.  Right  hand  side  of  the  linear  system  to  be 


solved . 


FI 


output.  The  solution  is  returned  in  the  array  FI  . 


N 


input.  Total  number  of  variables. 


NY 


EPS 


input.  Number  of  differential  and  nonlinear  variables, 
input.  Error  constant  from  LDASUB,  •  RMSEPS  . 


YMAX 


input.  Array  of  maximum  magnitudes  of  y^  up  to  the 
current  time  (or  one  if  maximum  magnitude  is  less  than  one). 


NEWPW  -  input.  Indicates  whether  a  new  J  matrix  has  been 


computed  since  the  last  entry  to  NUITSL. 


-  1,  indicates  this  is  a  new  J  matrix.  If  any 


preprocessing,  such  as  LU  decomposition  is  to  be  done, 


the  preprocessing  should  be  done  and  NEWPW  set  to  zero. 


=  0,  indicates  the  J  matrix  is  the  same  as  on  the  previous 


entry  to  NUITSL. 
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KRET  -  output.  Return  indicator. 

=  0  ,  normal  return 

=  1  ,  error  return,  solution  of  equations  not  obtained. 

Note  that  the  parameters  EPS  and  YMAX  are  useful  if  an  iterative 

method  is  used  for  solutuon  of  the  equations.  Because  the  solution 

represents  corrections  to  the  predicted  value,  and  corrections  to  that, 

the  solution  is  small  compared  to  the  dependent  variable  values.  Hence, 

compared  to  the  YMAX  array,  the  error  tolerance  can  be  fairly  large. 

The  following  convergence  criteria  have  been  used,  with  great  success. 

th 

Let  6u^  denote  the  i —  component  of  the  difference  between  successive 

th 

iterates,  with  u^  being  the  i —  component  of  the  current  iterate. 

Then  the  iteration  is  considered  to  have  converged  whenever 


Condition  (i)  requires  convergence  to  2  digits  more  accuracy  than  the 
user  has  asked  for  in  the  solution  of  the  system  (2.1),  relative  to 
YMAX.  Condition  (ii)  requires  the  same  relative  accuracy  in  u^  as  is 
asked  for  by  the  user  in  the  solution  of  the  system  (2.1),  unless  the 
solution  is  smaller  than  c  ,  in  which  case  the  change  is  compared  to 
c  rather  than  |u^|  .  This  avoids  difficulty  if  u^  is  close  to  zero. 
The  e  above  is  EPS  =  v^*RMSEPS,  where  M  and  RMSEPS  are  inputs  to 
SDESOL.  Two  versions  of  NUITSL  incorporating  iteration  and  this  conver¬ 
gence  test  are  given  in  the  appendix. 
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4.2.4  Subroutine  DERVAL 


This  subroutine  solves  for,  or  otherwise  supplies  the  initial  values 

of  the  derivatives,  and  possibly  other  variables.  In  some  instances 

it  may  need  to  be  supplied  by  the  user.  The  standard  version  of 

DERVAL  given  in  the  appendix  uses  Newton’s  method  to  solve  the  first 

m  (=NY)  of  the  equations  (2.1)  for  y(^Q)  *  assuming  values  for  y(t^) 

8F 

and  V(t^)  have  been  supplied.  To  accomplish  this,  the  matrix  — 

-20  ^ 

is  needed,  and  this  is  obtained  by  calling  JACMAT  with  h  =  16  , 

A2  =  -1  ,  and  N  =  NY  ,  implying  NL  =  0  for  this  call.  Special  care 


must  be  taken  if  in  fact  NL  is  not  zero  to  assure  that  the  matrix  is 

20  9F 

computed  and  stored  properly.  The  matrix  returned  is  then  16  —  . 

ay 


A  call  to  DIFFUN  yields  the  function  values 

F(y(to)  ,  yC^g)  9  ^q)  where  yC^g)  is  the  current  iterate. 

20 

Multiplication  of  the  function  values  by  16  and  a  call  to  NUITSL 

(again  with  N  =  NY)  gives  the  Newton  iterate.  Of  course,  the  same 

sort  of  special  care  as  necessary  in  JACMAT  is  necessary  in  NUITSL. 

3F 

Obviously  the  above  scheme  cannot  work  if  —  is  singular,  such 

9y 

as  it  would  be  if  one  of  the  equations  is  algebraic.  In  this  instance 
the  user  must  either  devise  his  own  version  of  DERVAL,  or  supply  the 
values  along  with  a  dummy  version  of  DERVAL.  In  an  extreme  case  the 
user  may  simply  set  initial  derivatives  to  zero.  This  will  provide  a 
poor  predicted  value  on  the  first  step,  and  will  force  an  artificially 
small  times tep  for  the  first  two  steps.  However,  the  overall  penalty 
is  generally  small,  as  appropriate  (corrected)  values  are  computed  at 
the  first  step,  and  after  two  steps  the  program  quickly  increases  the 
times tep. 
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The  calling  sequence  for  this  subroutine  is 
CALL  DERVAL  (Y,  YL,  T,  N,  NY,  DY,  KERET)  ,  where  the  parameters  are 
defined  as  follows. 

Y  -  input  and  output.  Same  as  in  SDESOL.  On  entry  Y(l,i) 

contains  the  initial  values  of  the  variables  .  On 

return,  the  values  of  the  derivatives  are  stored  in  Y(2  i) 
YL  -  input.  Same  as  in  SDESOL.  This  array  contains  the  initial 
values  of  the  linear  variables. 

T  -  input,  initial  time 

N  -  input.  Total  number  of  variables. 

NY  -  input.  Number  of  differential  and  nonlinear  variables. 

W  -  The  scratch  array  from  SDESOL,  can  be  used  in  any  way 

needed  by  this  subroutine. 

KERET  -  output.  Return  indicator 

=  0  normal  return 

=  1  error  return,  initial  values  were  not  obtained. 
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Appendix  1:  Program  Listings 


The  following  are  listings  of  the  basic  subroutine  package  and 
supporting  subroutines  which  are  of  general  use.  For  simple  problems 
the  user  only  needs  to  supply  a  calling  program  and  a  subroutine,  DIFFUN, 
to  evaluate  the  equations.  Use  of  the  NUITSL  routine  in  computer  facilities 
which  do  not  subscribe  to  the  IMSL  package  will  necessitate  modifications 
to  replace  LUDATF  with  another  LU  decomposition  routine,  and  LUELMF 
with  another  forward  and  backward  substitution  routine. 
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ooooooooooooooooooooooooooooooooooooooooooooooonooooooooooooooooooooooooo 


SLBROUTINE  SOESOL  < Y , YL iT ,T END, NY ,NL , M, JS K F  ,  M4X0E R , I PRT , H , HMI N 
IHRAX.RMSEPSfk) 


CALL 


SLEPQUTINE  SCESOL  !S  A  DRIVER  ROUTINE  FOR  SUBROUTINE  LDASU6, 
ITS  PURPOSE  IS  TO  SET  UP  THE  NECESSARY  REFERENCES  TO  A  LARGE 
BLOCK  OF  AUXILLARY  STORAGE,  AND  OBTAIN  INITIAL  VALUES  OF 
DERIVATIVES. 

THE  CALLING  SEQUENCE  FOR  SDESOL  IS 


WHERE  TFE  PARAMETERS  ARE  DEFINED  AS  FOLLOWS. 

Y  -  ARRAY  dimensioned  (7,NY>.  THIS  ARRAY  CONTAINS  THE 

DEPENDENT  VARIABLES  AND  THPIR  SCALED  DERIVATIVES. 
Y<J+l,n  CONTAINS  THE  J-TH  DERIVATIVE  OF  THE  I-TH  VA 
lABLE  TIMES  H*+J/ J-FACTGRI AL ,  WHERE  H  IS  THE  CURRENT 
STEPSIZE.  ON  FIRST  ENTRY  THE  CALLER  SUPPLIES  THE 
INITIAL  VALUES  OF  EACH  VARIABLE  IN  Y(1,I).  ON  SUB¬ 
SEQUENT  ENTRIES  IT  IS  ASSUMED  THE  ARRAY  HAS  NOT 
BEEN  CHANGED.  TO  INTERPOLATE  TC  NON-MESH  POINTS, 
THESE  VALUES  CAN  BE  USED  AS  FCLLOWS.  IF  H  IS  THE 
CURRENT  STEPSIZE  AND  VALUES  AT  TIME  T+E  ARE 
NEEDED,  LET  S  =  E /H  AND  THEN 


I-TH  VARIABLE  AT  T+E  IS 


JS 

SUM  Y( J  +  l  ,I  )*S+*J 
J  =  0 


YL 

T 

TEND 

NY 

NL 

M 

JSKF 


The  value  of  JS  IS  OBTAINED  IN  THE  CALLING  PPCGRAN 
BY  JS  =  IAbS( JSKF/101 

ARRAY  CF  NL  VARIABLES  WHICH  APPEAR  LINEARLY. 

CURRENT  VALUE  OF  THE  INCEPENCFNT  VARIABLE  (Tlvcj 
END  TIME 

NUMBER  OF  DIFFERENTIAL  EQUATICNS  AND  NONLINEAR 
VARIABLES. 

NUMBER  OF  LINEAR  VARIABLES 

NUMBER  OF  VARIABLES  INCLUDED  IN  THE  ERROR  TEST 
AN  INDICATOR  USED  8CTH  ON  INPUT  AND  OuTPL'^ 

CN  INPUT,  JSKF  =  -I  INDICATES  A  RESTART  CALL  TO 
SDESCL.  JSKF  *  0  INDICATES  AN'  INITIAL  CALL  TO 
SDESCL.  JSKF  >  0  INDICATES  A  CONTINUATICN  OF  THE 
PREVIOUS  CALL  TO  SDESOL.  JSKF  <  -I  MAY  HAVE  RESULTS 
FROM  THE  USER  NEGLECTING  TO  TEST  FOR  ERRCR  RETURNS 
FROM  SDESCL.  BECAUSE  OF  THIS  POSSIBILITY,  JSKF  <  -1 
RESULTS  IN  TERMINATION  CF  THE  RUN  WITH  THE 
APPROPRIATE  COMMENT. 

ON  OUTPUT,  JSKF  CONSISTS  OF  TwO  DIGITS  AND  SIGN, 

♦OR  -  QP.  Q  IS  THE  ORDER  OF  THE  FORMULA  CURRENTLY 
BEING  USED.  ?  INDICATES  THE  TYPE  OF  PETLRN,  AS 


FOLLOWS. 
JSKF  >  0,  F 
JSKF  <  0  IS 
MEANINGS, 


=  I  IS  THE  NORMAL  RETURN 

AN  ERROR  RETURN,  HITH  THE  FOLLOWING 


P  =  I  ERRCR  TEST  FAILURE  FOR  H  >  HMN 

P  =  3  CORRECTOR  FAILED  TC  CONVERGE  FOR  H  >  HM! 

P  =  4  CORRECTOR  FAILED  TC  CONVERGE  FOR  FIRST 

ORDER  METHOD 

P  =  5  ERROR  RETURN  FROM  SUBROUTINE  >IIIITsl 

P  =  6  ERROR  RETURN  FROM  SUBROUTINE  OERVAL 

MAXDER  -  maximum  ORDER  DERIVATIVE  THAT  SHOULD  BE  USED  IN 
METHOD.  IT  MUST  BE  NO  GREATER  THAN  SIX. 

IPRT  -  INTERNAL  PR’NT  CONTROL  INDICATOR  FOR  LDASU3. 

IPRT  =  0  NO  PRINT 

IPRT  >  0  PRINT  COUNTERS,  STEPSIZE,  CURRENT  TI(' 
AND  VALUES  OF  DEPENDENT  VARIABLES  AT 
EACH  STEP. 

H  -  CURRENT  STEPSIZE.  AN  INITIAL  VALUE  MlST  BE  SUPPLIEl 

BUT  NEED  NOT  BE  THE  ONE  WHICH  MUST  BE  USED,  SINCE  T^ 
SUBROUTINE  WILL  CHOSE  A  SMALLER  ONE  IF  NECCESSAPY  TC 
KEEP  THE  ERROR  PER  STEP  SMALLER  THAN  THE  SPECIFIED 
VALUE.  IT  IS  BETTER  TO  UNDERESTIMATE  THE  INITIAL 
STEPSIZE  THAN  TO  OVERESTIMATE  IT.  THE  STEPSIZE  IS 
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SCE 

30 
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4C 
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50 

SDE 

60 

SDE 

70 
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eo 
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SDE 
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SDE 
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SDE 
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SCE 

220 

SCE 

230 

SOE 

34C 
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SDE 

260 

SCE 

270 

SDE 

280 

SDE 

29C 
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SDE 
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SDE 
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SDE 
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SCE 

440 

OSDE 

450 

SCE 

460 

SCE 

470 

SOE 

480 

SCE 

490 

SCE 

500 

SCE 

510 

SDE 

520 

SCE 

52C 

SDE 

540 

SDE 

550 

SDE 

560 

SCE 

570 

N'SOE 

580 

SDE 

590 

SDE 

600 

SDE 

610 

SCE 

620 

SDE 

630 

SCE 

640 

SCE 

650 

SDE 

660 

esoE 

670 

SDE 

680 

SDE 

69C 

SDE 

7C0 

ESDE 

710 

720 

SDE 

730 

SOE 

74C 

SCE 

750 
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NORI^ALLY  NOT  CHANGED  PY  THc  LStR.  SDE  760 

HMIN  -  MINIf^UM  STEPSIZE  ALLOWED  770 

H.MAX  -  MAXIMUM  STEPSIZE  ALLOWED  -DE  780 

RMSEPS  -  THE  ERROR  TEST  CONSTANT.  THE  PIjCT-m«=  AN- SQUAR E  OF  SDE  7S0 

THE  SINGLE  STEP  ERROR  ESTIMATES,  ER  (  I  )  ,  OTVID’=D  3Y  800 

YMAX(I)  =  ('MAXIMUM  TO  CURRENT  TIME  OF  Yfll  )  MUST  3E  SDE  81C 

LESS  THAN  EPS.  THE  STcPSIZE  AND/OR  THE  CRDEP  SDE  82C 

ARE  VARIED  TO  ACHIEVE  THIS.  SDE  830 

W  -  SCRATCH  STORAGE  ARRAY.  MUST  EE  AT  LEAST  13*N!Y  ♦  5»NLSDE  8^0 

LOCATIONS,  PLUS  THOSE  REQUIRED  FOR  STORAGE  OF  THE  SDE  850 

MATRIX  PW  (SEE  DESCRIPTION  CF  SUBROUTINE  JACMAT).  SDE  860 

THE  STORAGE  JF  PW  WILL  NORMALLY  REQUIRE  NO  MORE  THAN  SCE  870 

*  2^N  LOCATIONS,  AND  IF  COMPACT  STnpAGF  TECH-  SCE  830 

NIQUES  AP«=  USED,  CAN  BE  MUCH  FEWER.  SDE  890 

SCE  900 
-SDE  910 
:0t  920 

SDE  920 
SOE  9A0 
SCE  950 
SCE  960 
SDE  97C 
SDE  960 
SDE  990 

soe  1000 


DIMENSICN  Y(7a),  YL(1),  Wd) 

IF  (JSKF.GT.O)  GO  TO  IZO 
IF  (JSKF.LT.-1)  GC  TD  140 
N  =  NY  +  N'L 

IF  (JSKF.LT.O)  GO  TO  110 

IF  THIS  IS  THE  FIRST  ENTRY,  OBTAIN  VALUES  CF  THE  DERIVATIVES. 
CALL  DFRVAL  ( Y, YL , T ,N, NY . W , KRET R ) 

IF  (KRETF.NE.O)  GC  TC  130 


THE  ARRAY 
THE  ARRAY 
THE  ARRAY 
THE  ARRAY 
THE  ARRAY 
THE  ARRAY 
THE  ARRAY 
THE  MATRIX 

110  NSVL  =  T’^NY  +  l 
NYMAX  =  NSVL^NL 
NER  =  NYMAX+NY 
NESV  =  .\EP  +  NY 
NFl  =  NESV+NY 
NCY  =  NFl+N 
NPW  =  NDY+N 
120  JS  =  JSKF 


^AGE 

BLOCKS  IN 

The 

W  ARRAY. 

THIS 

NF 

ECS  T 

0 

BE  CON'E 

SCE 

SOE 

1010 

1020 

AND 

SAVE 

ON  RESTARTS. 

STARTS  AT 

LOCATION 

1 

IN 

THE 

w 

ARRAY 

SDE 

SDE 

SDE 

1C20 

1040 

1050 

YLSV 

STARTS 

AT 

LOCATION 

NSVL 

IN 

THE 

w 

R  AY 

see 

1060 

YMAX 

ST^^nTS 

AT 

LOCATION 

N  YMAX 

IN 

THE 

w 

ARRAY 

set 

1070 

ER 

STARTS 

AT 

Lie  ATIGN 

NER 

1  N 

THE 

W 

AR'^  AY 

SCE 

1C8C 

ESV 

STARTS 

AT 

LOCATI jN 

ESV 

IN 

THE 

W 

APP  AY 

see 

1090 

FI 

STARTS 

AT 

locatijn 

NFl 

IN¬ 

THE 

w 

A  P  0  A  Y 

SCE 

1100 

DY 

STARTv 

AT 

LOCATION 

NCY 

IN 

THE 

w 

AY 

see 

1110 

PW 

STAR  TS 

AT 

LOCATION 

NPW 

I  N 

the 

w 

ARP  AY 

S  C£ 

1120 

CALL  LDASUe  (Y,  YL  , T  ,  T  FND  *  N  ,  NY  ,  ^  •  JS  ,  KF  •  M  AX DER  ,  I  PRT  ,  H  ,  ! N  ♦  HM  A  X  , 

IRMSEPS,  W, W<NSVL ) , w(NYMAX) , W ( Nr R ) , W ( NES V ) , W < N FI ) , W ( NC Y  )  , W ( NP W ) ) 


CODE  JSKF  ON  RETURN  FROM  LOASUB 

JSKF  ^  ISIGN(JS=<‘10+I  ABS(KF  )tKF) 
RETURN 

130  JSKF  =  -6 
RETURN 

140  PPINT  1,  JSKF 
STCP 


1  FORMAT  ('OIT  IS  AN  ERROR  T-O  ENTER  SDES; 
1  »  RUN  HAS  BEEN  TERMINATED.') 

END 


•L  WITH  JSKF  =  ' ,  IlD// 


SDE 

SCE 

SCE 

SDE 

SDE 

SDE 

SDE 

SCE 

SDE 

^CE 

SDE 

SCE 

SOE 

SOE 

SDE 

SDE 

SCE 

SDE 

SCE 

SCE 

soe 

SCE 

SCE 

SDE 

SDE 


1130 

1140 

1150 

1160 

1170 

1180 

1190 

1200 

1210 

1220 

1230 

1240 

125C 

1260 

1270 

1280 

1290 

1200 

1210 

1320 

1230 

1240 

1350 

126C 

1270 


SLBROUTINE  LCASUb  ( Y , YL , T , TEND ,M , NY , M t JST ART , KFLA G , M A XDR , I PRT,Ht 
IHM IN, HM AX, RMSEPS, SAVE, YL3V ,YMAX, ER, cSV, FI, DY,PW) 

SIBRQUTIK£  LCASUb  IS  A  MODIFICATION  OF  SUBRGUTINF  DFASUB 

WHICH  IS  DUE  TO  R.  L.  BROWN  AND  C.  W.  GEAR.  D'^ASUH  IS  DOCUMENTED 

IN  THE  REPORT 

C OCUMENTATIfN  FOR  0FA5UB-- 
BY  R.  L.  EROWN  AND  C.  W.  GEAR 
REPORT  UILCDCS-R-73-575,  JUlY  1973 
UNIVERSITY  OF  ILLINOIS  AT  UR EANA -C HAMP A  I 3N 
UR3ANA,  ILLINOIS  61801 
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LDA 
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THIS  REPCRT  IS  AVAILABLE  THE  .NATIONAL  TECHNICAL  INFORMATION 

SERVICE  CF  THE  U.  S.  DEPARTMENT  OF  COMMERCE  UNDER  ACCESSION  NUMBE 
CC0-IA69-225. 

THE  MODIFICATION  HERE  IS  DOCUMENTED  IN  THE  REPORT 

A  PRCGRAM  FOR  THE  NUMERICAL  SOLUTION  OF  LARGE  SPARSE  SYSTEMS  : 

ALGEBRAIC  AND  IMPLICITLY  DEFINED  STIFF  DIFFERENTIAL  EQUATIONS 

BY  RICHARC  FRANKE 

REPORT  NPS53FE76051t  MAY  1976 

NAVAL  POSTGRADUATE  SCHOOL 

MONTEREY,  CALIFORNIA  93940 


THE  CALLING  SEQUENCE  FOR  LDASUE  IS 

CALL  LDASUB  (Y,YL,  T,TEND,l^)>i'y,M.J.STAPT,KFLAG,MAXOR,  IPRT,H,HMINt 
hMAX,RMSFPS  ,S  AVE  ,  YLS  V  t  YMA)C,  E  R  ,  E$  V  ,  F 1 ,  DY  ,  P  W  ) 

k^HERE  THE  PARAMETERS  ARE  DEFI.NED  AS  FOLLOvJS* 

Y  -  ARRAY  DIMENSIONED  (7, NY).  THIS  ARRAY  CONTAINS  THE 

DEPENDENT  VARIABLES  AMD  TJIEIR  SCALED  DERIVATIVES. 
Y(J+1,II  CONTAINS  THE  j- fH  DERIVATIVE  OF  THE  I-TH  VA 
IABlE  TIMES  H^*J/J-FAC^ORIALt  ’^HERE  K  IS  THE  CURRENT 
STEPSIZE.  ?N  FIRST  ENTRY  THE  CALLER  SUPPLIES  THE 
INITIAL  VALUES  OF  EACH  VARIABLE  IN  Y(1,I)  AwD  AN 
ESTIMATE  OF  THE  InIITIAL  VALUES  CF  THE  DERIVATIVES 
IM  Y<2,I).  ON  SUBSEQUENT  ENTRIES  IT  IS  ASSUMED  THAT 
THE  ARRAY  HAS  NOT  BEEN  CHANGED.  TO  INTERPOLATE  T^ 
NCN-MESH  POINTS,  THESE  VALUES  CA.N  BE  USED  AS  FOLLOWS 
IF  H  IS  THE  CURRENT  STEPSIZE  AND  VALUES  AT  TIME  T+E 
NEEDED,  LET  S  =  E /H  AND  THEN 


I-TH  VARIABLE  AT  T+E  IS 


NQ 

SUM  Y( J+l  ,I  )*S*«J 
J  =  0 


THE  VALUE  OF  NQ  IS 
BY  NQ  =  JSTART. 


OBTAINED  IM  THE  CALLING  PROGRAM 


YL 


T 

TEND 

N 

NY 

M 


<0 


ARRAY  OF  NL  =  N  -  NY  VARIABLES  WHICH  APPEAR  LINEARLY 
THE  USER  SUPPLIES  INITIAL  VALUES  FOR  THESE  VARIABLES 
CURRENT  VALUE  OF  THE  INDEPENDENT  VAQZABLF  (TIME) 

END  TIME 

TOTAL  NUMBER  OF  VARIABLES 

NUMBER  CF  DIFFERENTIAL  EQUATIONS  AND  NONLINEAR 
VARI ABLES. 

NUMBER  OF  VARIABLES  INCLUDED  IN  THE  ERROR  TEST. 

THIS  NUMBER  CAN  BE  NO  GREATER  THAN  NY.  IF  IT  IS 
GREATER  THAN  NV,  NY  VARIABLES  ARE  USED  IN  ERROR 

TFST. 

JSTART  -  INPUT  AND  OUTPUT  INDICATOR. 

ON  INPUT  JSTART  HAS  THE  FOLLOWING  MEANINGS. 

THIS  indicates  a  RE-START  FROM  A  PREVIOUS 
POINT  FOLLOWING  TERMINATION  OF  THF  PUN  CR 
SOLUTION  OF  ANOTHER  PRCBLFM  DURING  1  HE  SAMF 
RUN.  PARAMETERS  IN  THE  CALLING  SEQUENCE 
MUST  HAVP  been  preserved  FROM  THE  PREVIOUS 
USE,  PARTICULARLY  THE  ARRAYS 
SAVE,  YLSV,  Esv,  AND  PW. 

THESE  ARRAYS  MUST  BE  SAVED  AFTER 
TO  SUBROUTINE  LDASAV,  WHICH  ALSC 
NECESSARY  PARAMETERS  INTERNAL  TC 
INDICATES  AN  INITIAL  CALL  TO  LDASU6.  THE 
ROUTINE  INITIALI^FS  ITSELF,  SCALES  thF 
DERIVATIVES  IN  Y(2,I)  AND  THEN  PERFORMS  THE 
INTEGRATION  UNTIL  T  >  TEND. 

INDICATES  THE  SOLUTION  IS  TO  BE  CONTIiNUED. 
AFTER  THE  INITIAL  ENTRY  IT  lo  NEITHER 
DESIRABLE  NOR  NCCESSARY  TO  RE-ENTER  WITH 
JSTART  =  0,  SINCfc  THIS  RE- I NI T I AL ! ZES 
THE  CODE,  BcGMNFNG  WITH  A  FIRST  CRDER 
METHOD  AGAIN. 

ON  OUTPUT,  JSTART  IS  SET  TO  THE  VALUE  OF  NC,  THE 
ORDER  OF  THE  FORMULA  CURRENTLY  BEING  USED. 


^0 


>0 


A  CALL 

SAVES 

LDASUB* 


LDA 

120 

RLDA 

130 

LDA 

140 

LDA 

150 

LDA 

160 

FLCA 
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LDA 

180 

LCA 

190 

LDA 

200 

LDA 
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LCA 
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LDA 
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-LDA 

240 

LCA 
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LDA 
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LDA 

270 

LCA 
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LCA 
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LCA 

300 

1  DA 
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LDA 

320 

LCA 
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RLDA 

340 

LOA 
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LCA 
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LCA 
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LCA 

380 

LCA 
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LCA 

400 

.LCA 

410 

LDA 
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LCA 
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LDA 

440 

LDA 
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LCA 

460 

LDA 

470 

LDA 

480 

LDA 

490 

LDA 
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LCA 

510 

.LCA 

520 

.LDA 

5  30 

LOA 

540 

LCA 

550 

LDA 

560 

LCA 

570 

LOA 

580 

LDA 

590 

LOA 

600 

LOA 

610 

LCA 

620 

LCA 

630 

LCA 

640 

LCA 

650 

LDA 

660 

LCA 

670 

LCA 

680 

LCA 

690 

LCA 

700 

LOA 

710 

LDA 

720 

LCA 

730 

LCA 

740 

LCA 

750 

LCA 

760 

LDA 

770 

LCA 

780 

LOA 

7S0 

LOA 

800 

LDA 

810 

LDA 

82C 

LCA 

8  30 

LDA 

840 

LDA 

850 

LCA 

860 
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KFLAG 


MAX?'5 

IP^T 


THE  CQNiPLETIGN 
MEANINGS 
+  1 
-1 
-3 
-4 


:00c  TNDICATU^,  WITH  THE  POLLUTING 


THE  INTEGRATION  WAS  5LCCESSFUL 
cRROR  TEST  FAILURE  FCR  H  >  HMN 
CORRECTOR  FAILED  TO  CCNVEPGE  FO-  H  > 

CORRECTOR  FAIlFC  to  CONVERGE  FQ'^  FIRST 
ORDER  METHOD 

-5  ERROR  RETURN  F‘=‘CM  S'JBkCUTIME  NUITSL 
MAXIMUM  ORDER  DERIVATIVE  THA^  SHOUuD  BE  LSED  IN  THE 
method.  it  must  be  mo  greater  than  Six.  if  IT  TS 
GRFATF.R  THAN  SIX,  THE  MAXIMUM  CRDEP  U'ED  NIlL  BE  SIX.LDA 


IC  A 
I.  CA 
LCA 
LDA 
LCA 
lCA 
LCA 
LCA 
LCA 
LCA 


HMIN 
HMAX 
RMS  EPS 


SAVE 

YLSV 

YMAX 


INTERNAL  PRINT  CONTROL  INPTCATCr 
^0  NO  PRINT 

>  0  PRINT  COUNTERS,  '7EFSI2E,  CURRENT  -^-TME 
AND  VALUES  OF  OE*'END£NT  VARIABLES  AT 

EACH  STEP.  _ 

CURRENT  STEP5IZE.  AN  iNITiAL  VALUE  MUST  BE  SUPPLIED  LDA 
BUT  NEED  MO-^  BE  THE  1\S  wH  CH  ^«ILL  6^  USED.  SIi^CE  ThELCA 
SUBROUTINE  WILL  CHCCSE  A  SMALLER  ONE  :F  NPCCESSARY  TCL DA 
KEEP  the  error  per  STEP  SMALLER  THAN  THE  SPECIFIEC  LCA 
value.  it  TS  eETTEP  TO  UNDERF ST  I  MAT E  THE  INITIAL 
S-^EPSIZE  than  TC  OVFkE:iTiMATE  IT.  THE  STEPSIZE  IS 
NORMALLY  NOT  CHANGED  BY  THE  USER. 

MINIMUM  STEPSIZE  ALLOWED 
MAXIMUM  STEPSIZE  ALLOWED 

THE  ERROR  TEST  CONSTANT.  THE  POOT-ME AN- SQUA ^ E  OF 
THE  SINGLE  STEF  ERROR  ESTIMATES.  £R (  I  )  ,  DIVIDED  EY 
YMAX(I)  ==  <MAXIMUM  TO  CURR'=NT  liME  OF  Y(T))  MUST  PE 
LESS  THAN  RMSEPS.  THE  STEOSIZF  AND/OR  CHDEF  APE 
VARI EO  TO  ACHIEVE  THIS. 

AN  ARRAY  0^  LENGTH  AT  LEAST  7*NY 
AN  4RRAY  CF  LENGTH  AT  LEAS"^  NL 

A  VECTOR  OF  LFNGTH  NY  WHICH  CONTAINS  THE  MAXTMyf/ 


LCA 

LCA 

LDA 

LCA 

LDA 


LCA 

LCA 

LDA 

LCA 

LDA 

LCA 

LDA 

LCA 

LCA 

LDA 

LCA 

LCA 

LDA 


OF  EACH  Y  SE*=N  SO  FAR.  CN  THE  FIRST  CALL,  THESF  WiLLLCA 


ER 

FSV 

PI 

DY 

PW 


BE  INITIALIZED  AS  YMAX(I)  =  M ^ X ( 1 , 1 Y (  1 ,  I  )  I) 

A  VECTOR  OF  LENGTH  MY 

A  VECTOR  OF  LENGTH  NY 

A  VECTOR  V  LFNGTH  N  =  NY 

A  VECTOR  jP  LENGTH  N  =  NY 

AN  ARRAY  IN  WHICH  THE  J 

IN  SUBROUTINE  JACMAT  WILL  BE  STORED.  SIZE  WHICH 
MUST  8£  ALLOWED  IS  DETERMINED  BY  THF  STCRAGE  TECH¬ 
NIQUE  USED  FOR  IT,  BUT  NOFMALIY  WON'T  BE  MORE  THAN 
^  2«M  LOCATIONS,  THE  LA'^’^rp  a^S'N  BEING  REQUIRED 
BY  THE  LINEAR  EQUATION  SOLVER. 


^  ML 
4-  NL 

MA'^PIX  CC-^PUTEC 


LCA 
LCA 
LDA 
LDA 
LCA 
LDA 
LCA 
LCA 
LCA 
LCA 
LCA 
LDA 

- 

OIMENSTZN  Y(7,l)»  YL<1),  5AVF(7,l),  YMAX(l).  ERd),  YLSV(l),  FKDLOA 

1.  PERT<A,3),  CCJF(21),  ESV  (  1  )  ,  DYd),  PW<1),  SAV(l).  A(29i  LDA 

EQUIVALENCE  (A(8),BND),  (A(9),BR),  <A(10),r),  ( A  <  1 1)  ,  ED  W  si )  ,  LCA 

l(A(12)f  ENQl),  (A(l3),FNQ2U  <  A  (  14  ) ,  ENP3  )  ,  (A<15),EPS),  (  A  (  16  )  ,  EUP  )  LD  A 

2,  < A< 17) ,HNtW ) ,  (A<18),PEPSH),  ( A ( 19 ) , I DOUB )  ,  ( A ( 20  )  , I  WE VAL )  ,  LCA 

3<A<21)fK),  (A(22) tLCOPYL) ,  ( A ( 2 3 ) , LCOP V Y ) ,  C A ( 24 ) , MA X CL R ) ,  LCA 

4<A(25),M1),  (A(26),NL),  (A(27),Ng),  (A<28),KSI,  (A(29),Nvy)  LCA 

LCA 

- LDfl 

LDA 

TrSTING  AND  LCA 

Dions.  LDA 

LCA 

- LCA 

»,l.,l.,.23,  LCA 

LCA 

- ICA 

LDA 

THE  STIFFLY  LCA 

MACHINE  LCA 

LCA 
LDA 
LCA 
LCA 
LDA 
LCt 
LDA 
lD  A 


TFE  COEFFICIENTS  IN  THE  PERT  ARRAY  ARE  USED  FjR  tRPCR 
CHANGING  STFPSIZE  AND  NEED  TO  BE  ACCURATt.  TC  CNLY  A  F 

CAT  A  PERT/4.,9.,16.,25.,36.,49.,9.,16.,25.,:-c.,49.,64 
12. 78d9E-2, 1.7C569E-3,6.83929E-5/ 

THE  ENTRIES  IN  THE  CCF  ARRAY  ARF  THE  COEFFICIENTS  FCR 
STABLE  METHCDS  USED  IN  THIS  PROGRAM  AND  ARE  TO  5c  THE 
PRECISION  EQUIVALENTS  OF  THE  FOLLOWING  CGNSIANTS. 

-1 

-3/2  ,  -1/2 
-11/6  ,  -1  ,  -1/6 
-25/12  ,  -35/24  ,  -5/12  ,  -1/24 
-137/60  ,  -15/8  ,  -17/2‘t  ,  -1/8  ,  -1/120 
-147/60,  -2C3/90  ,  -49/48  ,  -35/144  ,  -7/240  ,  -1/720 


no 

eao 

ESO 
900 
910 
S2G 
930 
S40 
9  50 
960 
970 
9  80 
990 
ICOO 
1010 
1020 
1033 
1040 
1C5G 
1060 
1C70 
1080 
1090 
1100 
1110 
1120 
1130 
1140 
1150 
1  160 
1170 
1180 
1  190 
1200 
1210 
1220 
1230 
1240 
1250 
1260 
1270 
1280 
129C 
1300 
1210 
1320 
1330 
1340 
1350 
1360 
13  70 
1580 
1390 
1400 
1410 
1420 
143C 
1440 
1450 
1460 
1470 
1480 
1490 
1500 
1510 
1520 
1530 
1540 
1550 
1560 
1570 
1580 
1590 
)6C0 
1610 


25 
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C  LDi  1620 

Q  - LDA  1620 

CiTA  C0F/-1.  ,-l  .  5, -.5 1  - 1.  8  33533,  - 1., 1666  667,-2. 0  833  23  » -1.458333,1-0  A  1640 

1- . 41666 6 7, -.04166 6  67, -2. 28  3  3 33, -1.8 75  ,- . 708233 3, - . 12 5  ,-.C0S 23 323 3 , LCA  16  50 

2-  2. 45, -2. 2  5  55  56, -1.0  206  33, -.243  0556, -.02  9  1666  7, -.00 13  83  8  39/  LD.4  16  60 

IF  (JSTA'STI  100,110,150  LDA  1670 

- LOA  1680 

IF  THIS  15  A  RESTART  ENTRY,  RESTORE  Y  A.'JO  VL  FROM  THE  S4VF  ANC  LCA  1690 

YLSV  ARRAYS,  ViHEP  E  THEY  HERE  SAVED  BY  A  PREVIOUS  CALL  TC  LOASAV.  LCA  1700 

- LOA  1710 

100  CALL  COPYZ  ( Y ,S A VE , LOOP YY )  LCA  1720 

CALL  COPYZ  ( YL, YLSV,LCOPYL )  LCA  172C 

GC  TO  150  LDA  1740 

- LOA  1750 

IF  THIS  IS  THE  FIRST  CALL,  INITIALIZE  YW.AX,  SCALE  CERIVATlVcS,  ANDLDA  1760 

INITIALIZE  INDICATORS  AND  SET  ORDER  TO  ONE.  LOA  1770 

FOR  OOUfiLE  PRECIS  ION,  SET  LCOPYL  =  14+NY  AND  LC3PYL  =  2*NL  IF  LOA  1780 

SLBRQUTINE  CCPYZ  IS  IK  SINGLE  PRECISION.  LCA  179C 

- LDA  1800 

no  KL  -  N-KY  LOA  1610 

LCOPYY  =  7*KY  LCA  1820 

LCOPYL  =  NL  LCA  1830 

R1  =  PIKC(R,KY)  LOA  1840 

EPS  =  SCRT  (  FLOAT!  HI  n*RKSEPS  LOA  1850 

MAXCEP  =  HIK0(MAXCR,6)  LCA  I860 

IF  (IPPT.LE.C)  GO  TO  120  LOA  1670 

PRINT  3,  N,KL,RMSEPS,TEND,H  LCA  1880 

PRINT  4  LOA  1890 

120  KS  =  0  LOA  1900 

NW  =  0  LOA  1910 

LCA  1920 

CO  130  J-1,KY  LOA  1930 

YHAX(J)  -  4HAXlIl.,AbS(Y(l, J)  ))  LOA  1940 

120  Y(2,J)  =  Y(2,J)*H  LDA  1950 

LDA  1960 

KC  =  1  LOA  1970 

eR  =  1.  LOA  1980 

ASSIGN  190  TC  IRET  LCA  1990 

- LOA  2000 

SET  COEFFICIENTS  FOR  THE  ORDER  CURRENTLY  BEING  USED.  LCA  2010 

E  IS  A  TEST  FOR  ERRORS  OF  THE  CURRENT  ORDER  NO  LCA  2020 

tUP  IS  TC  TEST  FOR  INCREASING  THE  ORDER,  ECWN  FOR  DECREASING  THE  LCA  2030 

CRCER.  LCA  2040 

- LDA  2050 

140  K  =  NCT  lKQ-1 )/2  LOA  2060 

CALL  COPYZ  (A(2),C0F(Ktn,N0)  LDA  2070 

K  =  NC+1  LCA  2080 

IC0U8  =  NC  LCA  2090 

ENCl  -  .5/NC  LCA  2100 

EKC2  *  .5/K  LCA  2110 

EKC3  =  .5/(KC+2)  LCA  2120 

PEPSH  =  EPS*»2  LOA  2130 

E  =  PERT (KC,1 IfPFPSH  LCA  2140 

cUP  =  PPRT(KC,2»*PEPSH  LCA  2150 

EOWN  =  PERT(KC,3)4PEPSH  LDA  2160 

BNC  =•  (  PPS*ENQ3)T*2  LCA  2170 

IREVAL  =  1  LOA  2180 

GC  TO  IRET,  (190,200,490,570)  LCA  2190 

150  IF  (H.FC.HNEW)  GO  TO  190  LCA  2200 

- LCA  2210 

IF  CALLER  HAS  CHANGED  H,  RESCALE  DERIVATIVES  TO  REFLECT  THAT  HNFW  LCA  2220 

HAS  USED  CN  THE  LAST  CALL.  LCA  2230 

- LDA  2240 

R  =  H/HKEH  LDA  2250 

assign  190  TO  IRET  LCA  2260 

GC  TO  61C  LCA  2270 

- LCA  2280 

SET  JS'^ART  TC  NQ,  THE  CURRENT  ORDER  OF  THE  HETHTL),  eSF3P=  EXIT,  LDA  2290 

AND  SAVE  THE  CURRENT  STEPSIZE  IN  HNFW.  LCA  2300 

-  LDA  2310 

160  JSTART  =  NQ  LOA  2220 

HNEW  =  H  LCA  2330 

RETURN  LCA  2340 

170  KS  =  NS+1  LCA  2250 

IF  (IPRT.LE.O)  GO  TO  180  I CA  2260 
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2CC 


PRINT  DATA  IF  DESIRED  BY  USER 

LDA 

PRINT  1,  NS,NW,NQ,H.T,(Y(ljI ) ,I=1,NY) 

IF  (NL.GT.O)  PRINT  2,  (YL ( I  )  , i^l,NL) 

CCNTINL-F 

IF  (KFLAG. LT.O)  GG  TO  16n 

I F  ( T.GE.TcNC )  GG  TO  160 

LCA 

LDA 

LCA 

LCA 

LDA 

TAKE  ANCTHtP  STEP  IF  T  <  TEND 

LCA 

JSTART  =  1 

ICA 

SAVE  DATA  FCP  TRIAL  WITH  A  SMALLER  TIMESTEP  IF  THIS 

STEP  FAILS 

LCA 

CALL  COPYZ  (SAVE, Y,LCGPYY) 

CALL  COPYZ  (YLSV, YL,LCaPYL ) 

RACUM  =  1. 

KFLAG  =  1 

HCLD  =  H 

NCOLD  =  NQ 

TCLD  =  T 

T  =  T+H 

HINV  =  l./H 

LDA 
LCA 
LCA 
LDA 
LCA 
LCA 
LDA 
LCA 
LCA 
—  1  r  A 

CCiViPUTP  FSeClCftO  VALUES  PY  SFFeCTTVPLY  MULTIPLYING 
VECTCR  PY  PASCAL  TPIANGL5  MA“PIX 

DERIVATIVE 

LCA 

LOA 

DC  210  J=2,K 

J3  =  K+J-1 

DO  210  J1=J,K 

J2  =  J3-J1 

LDA 

LCA 

LCA 

LCA 

LCA 

LDA 

OC  210  1=1, NV 

210  Y4J2,I)  =  Y(J2,n+Y(  J2  +  1,I  ) 


DC  220  1  =  1  , KY 
220  FR  (I  )  =  C. 


DC  270  L=l,2 

C&lL  DJPFUN  (  Y,  YL  ,T,HlNV,t)Y) 
IF  (I^PViSL.LT.l)  GC  TO  230 


IDA 
LCA 
LCA 
LOA 
LCA 
LDA 
LCA 
LCA 

- LDA 

DC  JO  TC  7HPEE:  CCfi'^ECTCk  ITERATIONS.  CGNVERGcNCE  IS  :BTAIN'FC  WHENLDA 

CHANGES  APE  LESS  THAN  6ND  ^HICH  IS  DtP‘.=  N0ENT  ON  THE  EPkCF  “^EST  LCA 
CONSTANT,  THE  SJH  OF  COR'^ECTIONS  15  ACCUMULATED  IN  Efi(I).  IT  IS  LCA 
EQUAL  T:  THE  K-TH  CERIYATIVE  OF  Y  TIMES  H*=i^K  /  (  K-F  AC  TO  R  I  AL  *A  (  K  )  ) ,  LDi 
AND  THUS  IS  PPOPOPT  ION?L  TO  THr  ACTUAL  EFRONS  TO  THE  L0r,-5T  PGV^F^  LCA 
PF  H  PRESENT,  WHICH  IS  H^-^'K.  LCA 

- lda 

! 

LCA 
LCA 
LDA 
-LCA 
LCA 
LCA 

LCA 

L  C  A 

LCA 
-LCA 
lCA 
LCA 
LCA 
LCA 
LCA 
LCA 
LCA 
LCA 
LOA 
LCA 
LCA 
LDA 
LCA 
LOA 


IF  THERE  HAS  BEEN  A  CHANGE  CF  ORDER  OR  THERE  HAS  BEEN  TRjUPLE 
WITH  CONVERGENCE,  PW  IS  R  E -b  V  ALUATE’j  PRIOR  TO  STARTING  th- 
CCRPECTCR  ITERATTCN.  IWEVAL  IS  THEN  SET  TC  -1  AS  AN  INCICATGR 
THAT  IT  HAS  BEEN  CON‘=.  NEWPW  IS  SET  ^':NZER^.  Tj  iNClCATE  TO 
^LBROUTINE  KUITSL  THAT  A  NEW  ?W  HAS  BEEN  PPCVlD'^O. 


CALL  JACRAT  ( Y , Yl , T , Hi NV , A ( 2 ) , N , NY , E PS , DY , F 1 , PW ) 
KFLAG  =  1 
IWEVAL  =  -1 
NW  =  NW+1 
N  E  W  P  W  =  1 

230  CALL  NLITSL  (  PW  ,  DY  ,  FI ,  N,r4Y  ,  E  PS  t  YM  AX  ,NFwP  w  ,  Kpn  cr  ) 
IF  (KRRET.NE.O)  GC  TO  bOO 
IF  (NL.Lc.O ) • GO  TO  250 

CO  240  1^1, NL 

240  YL  (I)  =  YL  (  I  )-Fl(  I  +  NY ) 


25C  CONTINUE 
CEL  =  0. 


2  3  70 
2300 
23RC 
2400 
2410 
2420 
2420 
2^40 
2450 
2460 
2470 
2480 
24S0 
2500 
2510 
2520 
2530 
2540 
2550 
2560 
2570 
2580 
2590 
2  6  00 
2tlC 
2620 
2630 
2640 
2650 
2660 
2670 
2680 
269C 
2700 
2710 
2720 
2730 

2  7  4C 
27CC 
2760 
2770 
2?B0 
2790 
2800 
2810 
2820 
2830 
2840 
2850 
2860 
2870 
2680 
2690 
2900 
2910 
2920 
2930 
2940 
2950 
2960 
2970 
2980 
2990 
3C0G 
3010 
3020 
3C30 
5040 
3050 
2060 
3070 
3080 
3090 

3  100 
3110 
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260 


270 


DO  260  1=1, N\ 

Y(  1,1)  =  Y  1,  n-FK  I  ) 

Y(2,I  )  =  Y(2,I)>A(2)«F1( I ) 

ER(1)  =  ER(I)+F1( I) 

DEL  =  DEL  +  ( FI (I )/AMAXl( YMAX( I  ),ABS( Y( 1, I ) ) ) )^*2 
CONTINUE 

IF  (L.GE.2)  BR  =  AMAXl ( .9* BR , DE L/DEL 1 ) 

DELI  =  CEL 

IF  (AMIM(0EL,BR^DEL*2.).LE.BND)  GO  TO  330 
CONTINUE 


280 

290 

200 

310 


T  =  TOLC 

IF  (IWEVAL)  280,300,290 

IF  (H.LE.HMIN^l. 00001)  GO  TO  310 

RACUM  =  RACING. 25 

CONTINUE 

GC  TC  560 

KFLAG  =  -3 


RESTORE  Y  ANC  YL  AFTER  CONVERGENCE  FAILURE 


320 

330 

340 


call  COPYZ  (Y,S AVE,LC0PYY) 
CALL  COPYZ  (YL, YLSV,LC0PYL ) 
H  =  HOLC 
NC  =  NQCLD 
GO  TO  170 


THE  CORRECTCR  CONVERGED,  SO  NC^^  THE  ERROR  TEST  IS  MADE. 


LDA  3120 
LDA  3130 
LD&  3140 
LDA  3150 
LDA  3160 
LDA  3170 
LDA  3180 
LDA  3190 
LDA  3200 
LDA  3210 
LDA  3220 
LDA  3220 
LDA  3240 
-LDA  3250 

TFE  CCRRECTIQR  ITERATION  FAILED  TO  CONVERGE  IN  3  TRIES.  VARIOUS  L CA  3260 
POSSIBILITI ES  ARE  CHECKED  FOR.  IF  H  IS  ALREADY  HMIN  AND  PW  HAS  LCA  3270 
ALREADY  BEEN  RE-EVALUATED.  A  NO  CONVERGENCE  EXIT  IS  TAKEN.  LDA  3280 
CTHERWISE  TFE  MATRIX  PW  IS  RE-EVALUATED  AND/OR  (IN  THAT  ORDER)  THELDA  3290 
STEP  IS  REDUCED  TO  TRY  AND  GET  CONVERGENCE.  LDA  3300 

-LDA  3310 
LDA  2220 
LDA  3230 
LDA  3340 
LDA  3250 
LDA  3360 
LDA  3370 
LDA  3380 
-LDA  3590 
LDA  3400 
-LDA  3410 
LDA  3420 
LCA  3430 
LDA  3440 
LOA  3450 
LC.A  3460 
-LDA  3470 
LDA  3480 
-LDA  3490 
LDA  3500 
LDA  3510 
LCA  3520 
LCA  3530 
LOA  3540 
LDA  3550 
LOA  3560 
LOA  3570 
-LDA  3580 
LDA  3590 
LDA  3600 
LDA  3610 
LDA  3620 
LDA  3630 
LOA  3640 
-LDA  3650 
LDA  366C 
LDA  3670 
LDA  3680 
LDA  3690 
LDA  3700 
LDA  3710 
LDA  3720 
LDA  3730 
LDA  3740 
LCA  3750 
LDA  3760 
LDA  3770 
-LDA  3780 
LDA  3790 
LDA  3800 
LDA  3810 
LDA  3820 
-LDA  3830 
LDA  3840 
LDA  3850 
LDA  3860 


C  =  0. 

DC  340  1=1, M 

YH  =  AMAXl ( AES(Y( 1, I ) ) ,YMAX( I ) ) 
D  =  D*( EP( I )/YM)**2 

IWEVAL  =  0 

IF  (C.GT.E)  GO  TO  380 


TFE  ERRCR  TEST  IS  OKAY,  SO  THE  STEP  IS  ACCEPTED.  IF  lOOUR 
NCW  BECOMES  NEGATIVE,  A  TEST  IS  MADE  TO  SEE  IF  THE  STEP  SIZE 
CAN  BE  INCREASED  AT  THIS  ORDER  OR  ONE  HIGHER  CR  ONE  LCWER. 

THE  CHANGE  IS  MADE  ONLY  IE  THE  STEP  CAN  BE  INCREASED  BY  AT 
LEAST  101.  ICOUB  IS  SET  TO  NC  TO  PREVENT  FURTHER  TESTING 
FOR  A  while.  if  no  CHANGE  IS  MADE,  I DOUB  IS  SET  TO  9. 

IF  (K.LT.3)  GO  TO  360 

DC  350  J=2,K 

DO  350  1=1, NY 

350  Y(J,I)  =  Y(J,I)-^A(J)^ER(I) 

26C  KFLAG  =  1 

ICOUB  =  IDCLB-l 
IF  (IDOLB)  410,370,510 
370  CALL  COPYZ  (ESV,ER,M1) 

GC  TO  510 

THE  FRRCR  TEST  FAILED.  IF  JSTART  =  0,  THE  DcPIVATIVcS  IN  THE 
SAVE  ARRAY  ARE  UPDATED.  TESTS  ARE  THEN  MADE  ^0  FIX  THE  STEPSIZE 
AND  PERHAPS  REDUCE  THE  ORDER.  AFTER  RESTORING  AND  SCALING  THE 
Y  VARIABLES,  THE  STEP  IS  RETRIED. 


380  IF  ( JSTART. GT.O)  GO  TO  400 
DO  390  1=1, NY 
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390  S^VE(2tI  )  =  \(2,I  ) 


400 

410 

420 

4  30 

440 

45C 

460 

470 


KFLAG  ^  KFLAG-2 
IF  (H.LE.HMIN)  GO  TO  550 
T  =  TOLC 

IF  (KFLAG,LE.-5 )  GO  TO  530 
PR2  =  ( C/E) ^♦ENQ2*1.2 
L  *  0 

IF  (NQ.LE.  I  )  GO  TC  430 
C  =  0. 

CC  420  J=lfM 

Y^^  =  AMAXK  AeS(Y(  I,  J  )  )  ,  YMAX(.J  )  ) 

C  =  0+( Y<K,J)/YM)**2 

PRl  =  (C/ECUK)*’>6NQ1*1.3 
IF  (PR1.GE.PP2)  GO  TO  430 
PR2  *  PRl 
L  -1 

IF  (KFLAG. LT.C. OR. NQ. GE. MAXOER  )  GO  TO  450 
C  =  0 

DO  440  J=lfM 

YP  =  AMAXl  (Ae$(Y(  1,  J)  ),Y^*AX(  J)) 

C  =  0+(  (Eft(  J)-ESV(  J)  )/YM)«’’'2 

PRl  =  ( C/EUF )*«ENC3*1.4 
IF  (PR1.GE.FR2)  GO  TO  450 
PP2  =  PRl 
L  =  1 

R  =  l./ANAXl(PR2f  1.5-5) 

IF  (KFLAG. LT.O. OR. R. GE. 1. 1 ) 

ICOUB  =  9 
GO  TO  510 
NEhC  =  NC+L 
K  =  NEWC-H 

IF  (NEV^C.LE.NQ)  GO  TO  480 
R1  =  A( KEWQ  )/FLOAT(NEWQ) 


GO  TO  460 


CC  470  J=ltNY 
Y(K,J)  =  ER(J)*R1 


480  CONTINUE 


IF  THF  STEP  WAS  OKAY,  SCALE  THE  Y  VARIABLES  IN  ACCORDANCE 
WITH  THE  NEW  VALUE  OF  H.  IF  KFLAG  <  0,  HOWEVER,  USE  THE 
SAVED  VALUES  (IN  SAVE  AND  YLSV).  IN  EITHER  CASE,  IF  THE  QRJER 
HAS  CHANGED  IT  1$  NECESSARY  TO  FIX  CERTAIN  PARAMETERS  BY  CALLING 
THE  PROGRAM  SEGMENT  AT  STATEMENT  NUMBER  140. 

ICDUB  =  NC 

IF  (N5wC.5Q.NC)  GC  TO  490 
NC  -  NEWC 

ASSIGN  490  TC  IRET 
GO  T<^  140 

490  IF  (KFLAG. GT.O)  GO  TO  500 
RACUM  =  PACUM^R 
GC  TO  560 

50C  R  =  AMAX1(AMIN1(HMAX/H,R) ,HMIN/H) 

H  =  H*R 
UEVAL  =  1 
ASSIGN  510  TC  IRET 
GC  TO  610 

510  CC  520  1=1, M 

520  YMAX(I)  =  AHAXl (AeS( Y( 1,1 ) ) ,YMAX(I ) ) 

GC  TO  170 

THE  ERROR  TEST  HAS  NOW  FAILED  THREE  TIMES,  SC  THE  DERIVATIVES  ARE 
IN  BAD  SHAPE.  RETURN  TO  FIRST  ORDER  METHOD  AND  TRY  AGAIN.  fF 
CCURSEf  IF  NC  =  1  ALREADY,  THEN  THERE  IS  NC  HOPE  AND  WE  EXIT  WITH 
KFLAG  =  -4. 

530  IF  (NQ.EC.l)  GO  TO  540 


LCA 

LCA 

LDA 

LCA 

LDA 

LDA 

LCA 

LCA 

LDA 

LCA 

LDA 

LCA 

LCA 

LDA 

LCA 

LCA 

LDA 

LCA 

LDA 

LDA 

LCA 

LDA 

LCA 

LDA 

LCA 

LCA 

LDA 

LCA 

LCA 

LDA 

LCA 

LOA 

LDA 

LCA 

LDA 

LCA 

LOA 

LDA 

LCA 

LCA 

LCA 

LOA 

LDA 

-LDA 

LOA 

LDA 

LCA 

LCA 

LDA 

-LUA 

LCA 

LCA 

LCA 

LCA 

LDA 

LCA 

LDA 

LCA 

LCA 

LCA 

LCA 

LCA 

LDA 

LOA 

LCA 

LDA 

LDA 

LDA 

-LDA 

LCA 

LDA 

LDA 

LCA 

-LDA 

LDA 


3  8  70 
3880 
2890 
3900 
3910 
3920 
3920 
2940 
2950 
3960 
29  70 
3980 
5990 
4000 
4010 
4020 
4030 
4040 
4050 
4060 
4070 
4080 
4090 

4  100 
4110 
4120 
4130 
4  140 
4150 
4160 
4  170 
4180 
4190 
4200 
4  210 
4220 
4230 
4240 
4250 
4260 
4270 
4280 
4290 
4  300 
4  310 
4  3  20 
4230 
4340 
4350 
4360 
4270 
4380 
4390 
4400 
4410 
4420 
4430 
4440 
4  450 
4460 
4470 
4480 
4490 
4  500 
4  5  1C 
4520 
4530 
4540 
^550 
4560 
4570 
4580 
4  590 
4600 
4610 
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NC  =  1 
ICOUB  =  1 

ASSIGN  570  TC  IRET 
GC  TO  lAO 
540  NCOLD  =  1 
KFLAG  =  -4 
GO  TO  220 
550  KFLAG  =  -1 
GC  TC  170 


THIS  SeCTirN  RESTORES  THE  SAVED  VALUES  OF  Y  AND  YL  ,  SCALING  ^ 
Y  DERIVATIVES  AS  NECESSARY,  AND  THEN  REiURKS  TO  THE  PREDICTOR 


560  H  =  HOLC’frRACLM 

H  =  AMAXKHMN,  ANIN1(H,HMAX)  ) 
570  RACUM  =  h/hCLD 
R1  =  1. 

CO  580  J=2tK 
R1  =  Rl’i'RACLM 


58C 


DC  580  I=1,KY 

Y(  J,n  =  SAVE(J,I  )=^R1 


CG  5S0  1=1, NY 
590  Y(l,n  =  SAVF(1,I  ) 

CALL  COPYZ  (YL,YLSV,LCGPYL) 
IREVAL  =  1 
GC  TO  20C 
6CC  KFLAG  =  -5 
GC  TO  16C 


THIS  SECTION  SCALES  THE  Y  DERIVATIVES  BY  R**J 


6  10  R1  =  1. 

CC  620  J=2,K 
PI  =  R1«P 

CC  620  1=1, NY 
620  Y(  J,n  =  Y  (  J,  n*Ri 

GC  TO  IRET,  (190,510) 


V 

OF 

E 


THIS  SECTION  ALLOWS  FOR  RESTARTS  AFTER  -OLVING  ANOTHER  PROPLE 
HAVING  TeP^«I^ATEO  THE  CURRENT  COMPUTER  kUN.  SUBROUTINE  LDA  SA 
SAVES  THE  NECESSARY  VALUES  WHICH  ARE  INTERNAL  TO  LCASLB.  FOR 
CCU6LE  PRECISION.  WITH  COPYZ  IN  SINGLE  PRECISION.  THE  NUMBER 
LCCATIONS  TC  PE  SAVED  J^ND  RESTORED,  LCOPYS  AND  LlOPYR.  MUST  fc 
SET  TC  56. 

IT  IS  ASSUMED  THAT  IN  ADDITION  TO  THE  V4RIAPLES  IN  THE  ARRAY 
SAVED  BY  CALLING  LDASAV,  THE  USER  ALSO  SAVES  THE  ARRAYS  SAVE, 
YLSV,  YHAX,  ESV,  AND  PW. 

TC  RESTART  THE  USER  FIRST  CALLS  LDARST  TC  RESTORE  THE  VALUES 
BY  LDASAV,  THEN  RE-ENTERS  LDASUB  WITH  JSTAPT  <  C,  AND  WITH  TH 
CTHER  PARAMETERS  THE  SAME  AS  RETURNED  FROM  THE  LAST  ENTRY  Tj 
LCASUB,  PARTICULARLY  THOSE  ARRAYS  MENTIONED  ABOVE. 

ENTRY  LCASAY(SAV) 

LCOPYS  =  29 

CALL  COPYZ  (SAV, A , LCOPYS » 

CALL  COPYZ  (SAVE, Y,LCOPYY ) 

CALL  COPYZ  (YLSV, YL,LCOPYL ) 

RETURN 

ENTRY  LCARST(SAV) 

LCCPYR  =  29 

CALL  COPYZ  (A,SAV,LCOPYR) 

RETURN 


LDA 
LCA 
LDA 
LCA 
LDA 
LDA 
LDA 
I  CA 
LDA 
-LDA 
HE  LDA 
LTDPLDA 

- LCA 

LCA 

LDA 

LCA 

LDA 

LCA 

LCA 

LCA 

LDA 

LDA 

LDA 

LCA 

LDA 

LDA 

LCA 

LDA 

LLA 

LDA 

LDA 

LDA 

LDA 

- LDA 

LCA 

- LDA 

LCA 

LDA 

LDA 

LCA 

LDA 

LDA 

LDA 

LDA 

LCA 

- LDA 

ORLDA 
LDA 
LDA 
LCA 
LDA 
LDA 
LDA 
LCA 
IDA 
LDA 
SAVECLDA 
LDA 
LDA 
LDA 
-LCA 
LDA 
LDA 
LDA 
LDA 
LCA 
LCA 
LDA 
LDA 
LDA 
LDA 
LDA 
LDA 
-LDA 
LDA 


462C 
4620 
4640 
4650 
4660 
467C 
468C 
4690 
4700 
4710 
4720 
4  7  30 

4  7  40 
4750 
4760 
4770 
4780 
4790 
4800 
4810 
4820 
4830 
4840 
4850 
4860 
4870 
4880 
4690 
4900 
4910 
4920 
4930 
494C 
4950 
4960 
4970 
4980 
4990 
5000 
5010 
5020 
5030 
5G40 
5050 
5060 
5070 
5080 
5090 
5100 
5110 
5120 
5130 
5140 
5150 

5  160 
5170 
5180 
5190 
5200 
5210 
5220 
5230 
5  240 
5250 
526C 
5270 
5280 
5290 
5200 
5310 
5320 
5230 
5340 
5350 
5260 
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1  FCf5f^Ar  (2I5,I2»  lo2E10.2,7El4.6/(  32X,7E14.6)  i 

2  FORMAT  (:-2X  ,1F7E14.6) 

3  FORMAT  l*l  N  ='.I3,»  NL  =M3,*  FM^EPS  =*,lP69.2t» 

1  it9.2,'  F  =SE9.2//) 

4  FORMAT  (•  NS  NW  Q  H'»8X,*T 

PNC 


TuNC 


•t8X,»Yil,«)  AMC  VL(^)’//) 


LCA  3270 
LCA  5280 
LCA  5290 
LCA  5400 
LQA  5^10 
LCA  5^20 
LCA  5430 


SL8R0UTIN0  CCPYZ(StY,L) 
DIMENSION  S(l),Y(lJ 


This  SUBROUTINE  COPIES  THE  ARRAY  Y,  OP  LENGTH  L,  INTO  TPP  AnFAY 


IF (L.Lc.C  )RtTURN 
CG  100  J=ltL 
IOC  S(J)  =  Y(j} 
RETURN 
PNC 


COP 
COP 
---CC-P 
rcp 
S  COP 
CCF 
■--OOP 

COP 

COP 

COP 

CCF 


10 

20 

30 

40 

50 

60 

70 

80 

90 

LOO 

110 

120 


SLBROUTINE  CERVAL  ( Y t YL t T t N t N Y , U t KF RET ) 


DER 
DGR 
CEO 


C 

c 


THIS  SUBRCUTINF  CALCULATES  THE  INTIAL  VALUES  uH  THE  CFRTVAT*VPS 
IN  THE  GENERAL  CASE.  IT  IS  WRITTEN  SO  THAT  IT  SHOULD  WORK  IF  THP  CEP 
FIRST  NY  EOLATIONS  ALL  INVOLVE  DERIVATIVES.  1“  ATTEMPTS  TO  SOLVE  CEP 
THE  FIRST  NY  EQUATIONS  USING  NEWTON'S  METHOD,  BUT  SINCE  IT  TRIES  CEB 
TC  EVALUATE  CF/DY'  BY  CALLING  JACMAT  IN  SUCH  A  WAY  AS  TO  MAKL  the  OER 
CF/DY  TERM  INSIGNIFICANT,  IT  IS  POSSIBLE  THAT  IT  MAY  FAIL  FOP  tha-^OEP 
REASON.  IT  MAY  FAIL  FOR  OTHER  REASONS,  AS  WELL.  IF  IT  DOES  FAIL  CEP 
THE  USER  CAN  SUPPLY  HIS  uWN  VERSION  OF  DERVAL,  OP  MODIFY  TH.5  CEP 

ROUTINE  IN  SUITABLE  FASHION.  THIS  ROUTINE  ASSUMES  THAT  VALUES  CF  CEP 
THE  linear  VARIABLES  HAVE  BEEN  SUPPLIEC  PREVIOUSLY.  IF  THOSE  CFR 

MUST  BP  SCLVEC  FOR  SIMULTANEOUSLY  WITH  THE  CEPiVATIVtS,  THF  USER  DER 
VLST  SUPPLY  HIS  OWN  VERSION  OF  DERVAL.  DFP 

OPR 
DEF 
CEP 
CER 
OtR 
CER 
DER 
DEO 
CEP 
CER 


THE  CALLING  SEQUENCE  FOR  THIS  SUBROUTINE  IS 

CALL  D'^PVAL(Y,YL,T,N,N'r,W,KCRET  J 

WHERE  THE  PARAMETERS  ARE  DEFINED  AS  FOLLOWS 


YL 

T 

N 

NY 

W 

KERET 


SAME  AS  IN  LDASUB  AND  SOESOL.  Y(l,n  CONTAINS  THE 
INITIAL  VALUES  OF  THE  DEPENDENT  VARIABLES.  THF 
VALUES  OF  THE  DERIVATIVES  ARE  OETURNED  IN  Y(2,^). 

SAME  AS  IN  LDASUB  ANC  SOESOL.  THE  INITIAL  VALJCS  OF  CEP 
THE  LINEAR  VARIABLES  MUST  BE  SUPPLIED  TO  THIS  yrRSIClNuEF 
initial  TIME  DEP 

SAME  AS  IN  LDASUB,  TOTAL  NUMBER  OF  VARMELES  CER 

SAME  AS  IN  LDASUB,  NUMBER  OF  CIFFFRENTIAL  EQUATIONS  DEP 
AND  NONLINEAR  VARIABLES  DEP 

SCRATCH  ARRAY  W  FROM  THE  CALLING  SEQUENCE  OF  SDESCL.  CtP 
THIS  CAN  BE  USED  AS  NEEDED  IN  this  SUBROUTINE.  DEP 

RETURN  INDCATOP  DER 

=0  NORMAL  PETURN  DER 

=1  ERROR  PETURN  DER 

DER 

- CER 

DEB 
CcR 
DER 
DPR 
Dtp 
DER 
DEP 
PER 
DER 
DtR 
DEP 
DER 
DEP 
DcR 
DER 


DIMENSION  Y(7,l),  YL(1),  W(l) 

CC  100  1=1, NY 

W(2<‘N-H)  =  AMAXl(ABS<Y(l,n  1,1.) 
ICO  Y(3,I )  =  0. 

HINV  =  16. **20 
KERET  =  C 
EPS2  ^  NY/1.E8 
EPS  =  SCRT(EPS2) 

DC  140  IT  =  iaO 

DC  110  1=1, NY 
llO  Y(2,I )  =  Y(2,I)/HINV 


10 
20 
30 
40 
50 
60 
70 
80 
90 
100 
lie 
120 
12G 
140 
150 
160 
170 
180 
190 
200 
210 
220 
230 
240 
250 
260 
27C 
280 
290 
200 
310 
220 
3  30 
240 
250 
2  60 
270 
28C 
390 
^00 
410 
420 
420 
^40 
450 
^60 
470 
^80 
^90 
500 
510 
520 
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120 


130 


140 


150 

160 


170 


Ct^LL  DIFFUN  ( Y ,  YL  ,  T  ,  HI  NV  ,  W  ) 

Ci^LL  J£C^^AT  <Y,YL.T,HINV,-1.,NY,NY,EPS,W,W(N^1)  ,W(3’<‘N  +  1)  ) 
NEUPVi  -  1 

CC  120  I^lfNY 
U(  I)  =  V^d  )’*HINV 

CALL  NUIT5L  (W(3  +  N-H  ) »  W ,  W  (  N  + 1 ) ,  NYt  NY  ,  C  PS  ,  U<  2<^N+1)  »  NF  W  Fl^  ,  K  RE  7  ) 
IF  (KPPT.NF.O)  GO  TG  170 
ER  =  0. 

DC  130  1=^1, NY 

Y(3,n  Y(3,I)-W(N*I) 

ER  =  EP+(W(N+I)/AMAXl(ABS(Y(3fI  ))»1.)  )*<=2 

IF  (EP.LT.EPS2)  GO  TO  150 
CONTINUE 

GC  TO  17C 

CC  160  1=1, NY 
Y(2,I)  =  Y(3,I) 

RETURN 
KERET  =  1 
RETURN 
PNC 

SLBFOUT  INE  JACMAT  (  Y ,  YL  ,  T  ,  HI  NV ,  A2  ,N  ,  N  Y  , E P S  ,  DY ,  F  1 ,  P  ) 


SLBPCUTINE  JACMAT  IS  (USUALLY)  SUPPLIED  BY  THE  USER.  ITS  P 
IS  TO  EVALUATE  THE  J  MATRIX  NEEDED  WHEN  THE  CORRECTCR  EQU 
IS  solved  by  NEWTCN'S  METHOD.  THIS  VERSION  APPROXIMATES 
J  BY  NU^^ERICAL  CIFFERENCING  AND  USES  FULL  STORAGE  MOLE 
IN  AN  NXN  Matrix. 


URPCS5 

mTIGN 


JACMAI 


:alculates  the  matrix 


J  = 


DF  A2  DF 

CY  H  DY» 

THE  CALLING  SEQUENCE  FOR  THIS  SUBROUTINE  IS 
CALL  JACMAT ( Y.YLiTiHlNV,A2 ,EPS,N,NY,DY,F1,PW) 

where  the  parameters  are  defined  as  FOLLOWS. 

Y  -  SAME  AS  IN  LDASUB  AND  IN  SDESOL.  ON  INPUT  TO  FFIS 

SUBROUTINE  THE  ARRAY  CONTAINS  CURRENT  VALUES  or  THE 
DEPENDENT  VARIABLES  AND  THEIR  (SCALED)  DERIVATIVES. 

YL  -  SAME  AS  IN  LDASUB  AND  IN  SDESCl.  ON  INPUT  TO  THIS 

SUBROUTINE  THE  ARRAY  CONTAINS  CURRENT  VALUES  OF  THE 
LINEAR  VARIABLES. 

T  -  CURRENT  TIME 

HINV  -  1/H  ,  WHERE  H  IS  THE  CURRENT  STEPSIZE 

A2  -  A(2)  FROM  LDASUB. 

N  -  SAME  AS  IN  LDASUB,  TOTAL  NUMBER  OF  VARIABLES 

NY  -  SAME  AS  IN  LDASUB,  NUMBER  OF  DIFFERENTIAL  EOUhTIONS 

AND  NCNLINFAR  VARIABLES 

EPS  -  L2  ERROR  CONSTANT  USED  IN  LDASUB; 

OY  -  ARRAY  OF  FUNCTION  VALUES  A*  CURRENT  VALUES  OF  THE 

VARIABLES,  INPU'»*  TO  JACMAT- 

FI  -  SCRATCH  ARRAY  OF  N  LOCATIONS  WHICH  CAN  BE  USED  EY 

THIS  SUBROUTINE  IN  ANY  WAY  NEtCED. 

PW  -  J  MATRIX,  OR  APPROXIMATION,  CALCULATED  IN  JACMAT  AN 

RETURNED  TO  CALLING  PROGRAM.  THIS  MATRIX  IS  USED  IN 
SUBROUTINE  NUITSL  AND  STORAGE  MODE  MUST  AGREE  EETWEE 
THE  TWO  SUBROUTINES. 

DIMENSION  DY(1),  Y(7,l),  YL(1),  FlU),  °W(1) 


DER 

5  30 

DEP 

540 

DER 

550 

OER 

560 

DER 

570 

HER 

580 

DER 

5<50 

DEP 

600 

HER 

610 

DSR 

620 

DEP 

630 

DER 

640 

DcP 

650 

DER 

660 

DEP 

670 

DEP 

680 

DEP 

690 

DEP 

700 

DER 

710 

DER 

720 

DER 

730 

DEP 

740 

DER 

7  50 

DER 

760 

DER 

770 

DER 

780 

DEP 

790 

DER 

800 

JAC 

10 

-JAC 

2? 

JAC 

3C 

JAC 

40 

JAC 

50 

JAC 

60 

JAC 

70 

JAC 

80 

JAC 

90 

JAC 

100 

JAC 

no 

JAC 

120 

JAC 

130 

JAC 

140 

JAC 

150 

JAC 

160 

JAC 

170 

JAC 

180 

JAC 

190 

JAC 

200 

JAC 

210 

JAC 

220 

JAC 

230 

JAC 

240 

JAC 

250 

JAC 

260 

JAC 

270 

JAC 

280 

JAC 

290 

JAC 

300 

JAC 

310 

JAC 

320 

JAC 

330 

JAC 

340 

JAC 

350 

JAC 

360 

JAC 

5  70 

JAC 

330 

DJAC 

290 

JAC 

400 

NJAC 

410 

JAC 

420 

JAC 

430 

-JAC 

440 

JAC 

450 

32 


oonoooooooonoooooooooonooooonoooooo 


M  =  N-NV 
=  N^N 
C 

CC  100  1=1, NN 

100  Pt^d)  =  0. 

c 

c 

CC  12C  J=1,NY 
F  =  Y(  1,J J 
E  =  Y(2,J) 

=  EPS+^^^^X1(EPS,ABS(F),ABS(E)  ) 
Y(l,  J)  =  Y(  1,J)+R 
Y(2,  J)  =  Y(2,  J)-A2^R 
call  OIPFUN  (Y, YL  ,T,HINV,F1  ) 

C 

CC  110  1=1, N 

110  PVy(!  +  (  J-1)*N)  =  (  Fl(  I  )-0Y(  I  n/R 
C 

Y(2,  J)  =  F 
12C  Y(1,J)  =  F 
C 

TF  (M.  EC.O  )  GO  VO  150 
C 

CC  140  J  =  1,M 
F  ^  YL(J) 

R  =  FPS^AMAX1(EPS,A0S(F)  ) 

YL(J)  =  YL(JJ+R 

CALL  DIFFUN  ( Y, YL , T , HI N V , F 1 » 

C 

CC  130  1  =  1, N 

130  Pk(I+( J4NY-1)«N)  =  (Fl(I)-JY(I J )/R 
C 

140  YL(J)  =  F 
C 

150  CCNTIN'UE 
RETURN 
EMC 


J  AC 
JAC 
JAC 
lAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 


460 
470 
480 
490 
500 
510 
520 
530 
540 
550 
560 
570 
580 
590 
600 
6  10 
62C 
630 
64C 
650 
660 
670 
680 
690 
700 
710 
720 
730 
740 
750 
760 
770 
78C 
790 
6CC 
810 
820 


SUBROUTINE  NUITSL  (  PW,  OY ,  F 1 ,  N  ,NY  ,  EPS  ,  Yr'AX ,  NEWPW ,  KR  ET  ) 

THE  FUPFCSE  OF  THIS  SUBROUTINE  IS  TO  SOLVE  A 
LINEAR  SYSTEM  Of  EQUATIONS  FOP  THE  NEWTON  ITERATES  WHEN  THE 


NUI 

-NUI 

NUI 

NUI 


CORRECTuP  EQUATION  IS  BEING  SOLVED.  UPON  ENTRY  TO  THIS  SUB aC UT I N ENU I 


THE  SYSTEM  CF  EQUATION'S  TO  BE  SOLVED  IS 
J  IS  STORED  IN  PW  UPON  ENTRY 
W  I S  RETURNED  IN  FI 
-F  IS  STCRED  IN  DY  UPON  ENTRY 


J  W  = 


WHFRE 


NUI 

NUI 

NU! 

NUI 

NUI 


THIS  SUBROUTINE  IS  GENERALLY  SUPPLIED  BY  THE  USER,  ALTHOUGH  THERE  NUI 
ARE  SOME  STANCARO  FORMS  AVAILABLE.  FOR  EXAMPLE,  iHlS  VERSION 
ASSUMES  THAT  PW  IS  STORED  IN  FULL  STOPAC-E  MODE  IN  AN  NXN  MATRIX. 

IF  NEWPW  =  1,  AN  LU  DEC  GM  PO  5 !  T I CN  IS  DONE,  NEWPW  IS  SET  *^0  ZERO 
A^,C  FDRWARC  AND  BACKWARD  SUBSTITUTION  FOR  THE  SOLUTION  IS  DONE. 

IF  NEWPW  ==  0,  ONLY  FORWARD  AND  BACKWARD  SUBSTITuTTCN  FOR  THE 
SOLUTICN  IS  NECESSARY. 


NUI 
NLI 
NU! 
NUI 
NUI 
NUI 
NUI 

NCTE  THAT  THIS  VERSION  OF  NUITSL  REQUIRES  THAT  PW  HAVE  N^  +  2  +  2-r‘  NUI 
LCCATICNS  SINCE  2 ♦N  LOCATIONS  ARE  US^D  BY  THE  IMSL  LINEAR  EQLATIONNUI 
SCLVEP.  NLI 

NUI 

NCTE  that  THE  PARAMETERS  EPS  AND  YMAX  ARE  USEFUL  IF  aN  ITERmTIvE  NUI 
METHOD  IS  USED  TO  SOLVE  THE  SYSTEM  OF  EQUATIONS.  NUI 

NUI 

THE  CALLING  SEQUENCE  FOR  THIS  SL'BROU*riNE  IS  NUI 

NUI 

CALL  NUITSL (PW,DY  ,F 1 , N , NY , E P S , YMA X, NE WP W , KP ET )  NUl 

NUI 

WHERE  the  parameters  ARE  DEFINED  AS  FOLLOWS.  NLI 

NUI 

PW  -  THE  J  MATRIX  CALCULATED  IN  SUBRCUTINE  JACMAT  NUI 

DY  -  THE  RIGHT  HAND  SIDE  OF  THE  LINEAR  SYSTEM  TO  BE  SuLVEDNUl 

FI  -  THE  SOLUTION  IS  RETURNED  IN  THE  ARRAY  FI  NUI 

N  -  SAME  AS  IN  LDASUB,  TOTAL  NUMBER  OF  VARlABLtS  NUI 

NY  -  SAME  AS  IN  LDASUB,  NUMBER  CF  C I FFERE  N*^  I A  L  FQJATICNS  NUI 


10 

20 

30 

40 

50 

60 

70 

80 

90 

100 

110 

120 

130 

140 

150 

160 

170 

180 

190 

200 

210 

220 

230 

240 

250 

260 

270 

280 

290 

300 

310 

320 

^30 

340 

350 

360 
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EPS 

YMAX 

NEWPH 


KRET 


AND  NONLINEAR  VARIABLES 
L2  ERROR  CONSTANT  USED  IN  wDASUB 

MAXIMUM  VALUES  OF  Yd.I)  SEEN  UP  TO  THE  CURRENT  TIME 
INDICATES  WHETHER  A  NEW  J  MATRIX  HAS  BEEN  COMPUTED 
=1  INDICATES  A  NEW  J  MATRIX  HAS  BEEN  COMPUTED 
SINCE  THE  LAST  ENTRY  TO  NUITSL.  NEWPW 
SHOULD  BE  SET  TO  ZERO  IF  SOME  PREPROCESSING 
SUCH  AS  LU  DECOMPOSITION  MUST  BE  DONE  TN  A 
NEW  J  MATRIX. 

»0  INDICATES  THE  J  MATRIX  IS  THE  SAME  AS  WHEN 
NUITSL  WAS  LAST  ENTERED 
RETURN  INDICATOR 

=0  NORMAL  RETURN 

=^1  ERROR  RETURN.  SOLUTION  OF  EQUATIONS  CCULD 
NOT  BE  OBTAINED. 


CIMENSICN  PW<1),  DY(l)f  Fid),  YMAX(l) 

NL  =  N-NY 

IF  (NEWPW. EC. 0)  GO  TO  100 
NEWPW  =  0 
NN  *  N*  +  2«H 
NNN  *  NN+N 

CALL  LUCATF  (PW • P W , N.N . 0, D 1 , D2» PW (NN ) » PW ( NNN ) , F I , I ER ) 
IF  (  lER.EO.O)  GO  TO  io6 
KPET  »  1 
RETURN 

IOC  CALL  LUELMF  ( PW , D Y , PW ( NN) , N , N , F 1 ) 

KRET  =  0 

RETURN 

EKC 


NUI 

Nu: 

NUI 

NU' 

NU 

NUI 

Nu: 

NUI 

NUI 

NU 

NU_ 

NUI 

NU 

NU 

NU 

NU_ 

'NUI 

NU 

NU 

NU 

NU 

NU 

NU 

NU 

NU 

NUI 

NU 

NU 

NU 

NU 

NU 


SUBROUTINE  Cl FFUN ( Y , YLtT , HINV f DY ) 

ilBROUTINE  cTffuN  IS  SUPPLlio  BY  THE  USER.  ITS  PURPOSE  IS  TC 
EVALUATE  THE  FUNCTIONS  AT  CURRENT  VALUES  OF  THE  VARIABLES. 

THE  CALLING  SEQUENCE  FOR  THIS  SUBROUTINE  IS 

CALL  DIFFUN(V,YLiTtHINVf DY) 

WHERE  THE  PARAMETERS  ARE  DEFINED  AS  FOLLOWS. 

Y 

YL 
T 

HINV 
DY 


SAME  AS  IN  LDASUB  AND  SDESQL.  ON  INPUT  TO  THIS 
SUBROUTINE  THE  ARRAY  CONTAINS  CURRENT  VALUES  OF  THE 
DEPENDENT  VARIABLES  AND  THEIR  (SCALED)  DERIVATIVES. 
SAME  AS  IN  LDASUB  AND  SDESOL.  ON  INPUT  TO  THIS 
SUBROUTINE  THE  ARRAY  CONTAINS  CURRENT  VALUES  OF  THE 
LINEAR  VARIABLES. 

CURRENT  TIME 

1/H  .  WHERE  H  IS  THE  CURRENT  STEPSIZE 
RETURNED  ARRAY  OF  FUNCTION  VALUES. 


DIMENSICN  Y(7,1),YL(1),DY(1) 

DEFINE  YOLR  FUNCTION  HERE 

RETURN 

EKC 


270 
380 
290 
AOO 
410 
420 
430 
440 
450 
460 
470 
480 
490 
500 
510 
520 
530 
540 
550 
560 
570 
580 
590 
600 
610 
620 
630 
640 
650 
•  660 
670 
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Appendix  2:  Examples 


Example  1:  This  example  is  the  problem  proposed  by  Gear  [3].  The 
system  of  equations  is 


^  S 


+  (R-y^)^  +  I 

j=l  J  J 


=  1,  2,  3,  4 


1  r 

where  R  =  ^  L  y  • 
i=l  ^ 


1  ^  9 

s  =  i  I  (R-y.)^  , 

i=i 


and 


ys  +  y^  yg  +  y^  yg  =  0 


2yg  +  yg  -  yj^  +  -  1  -  e  =  0 


''l  -  ^^2  +  ^1  ^6  =  0 


+  V2  +  5y^  y2  =  0  . 


In  the  above  b--  =  b^^  =  b^«  =  b,,  =  447,501 
11  22  33  44 


'>12  ■  -  ‘’21  =  -‘’43  -  - 


‘>13  ■  -'>24  -  >>31  '  -‘’42  -  - 


‘'14  -  ‘>23  ‘41  -‘'32  -  • 


The  initial  conditions  are 

Vj  =  1  ,  1-1,  2.  3,  4  . 

yj  -  yo  -  1 

Vi  =  -  2  ,  V2  =  -  3  . 

9F 

Note  that  a  different  version  of  DERVAL  is  necessary  since  —  is  singular. 

9y 
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CIMCNST  :N  Y(7,6I,YL(2),W(1^)0),3!  (161 
/CfiT/G(16) 

DATA  G 1/4^7.501 ,-452.499,-47.499,-52. 501  , -452 . 499 , 447 .5 Ji , 52.501  , 

1  4  7.49  9,-4  7.499, 5  2.  501, 4*+?.  501,452.49  9,-52.501,47.499,452.4  99, 

2  447.501/ 

OAT  A  N,  NY,M  ,  ^,REPS  ,HMAX,  T,  TEND/8, 6 ,2,6, 1.6-2,  5.F.2,1.E-10, 

1  l.E-4,0., 1.F3/ 

CALL  EPPScT (207,256,-1,1 ) 

CALL  EPPSET  (208,256,-1,  1  ) 

CG  8  1^1,16 
8  G( I)  ^  GI (I  ) 

DC  10  1=1,4 
1C  Y( 1,1)  =  -1. 

Yd, 5)  =  1. 

Y(l|6)  =  1. 

YL( 1)  =  -2. 

YL(2)  =  -2. 

JSKF  =  0 

call  SOESGL  (Y,YL,T,T6N0,NY,NL,'^,  JSKF,6,  l,H,HWliN,HviAX,  ^EP'£,  W) 

PRINT  6, JSKF 

6  FCRMAT(  *0  JSKF  =*  ,14) 

STOP 

ENC 


SLBPGUTINE  CIFFUN  (  Y  ,  YL  ,  T  ,  h  I  .nIV  ,  DY  ) 

CC?MjN  /CAT/G(4,4) 

DATA  TGLC/-12.459/ 

0INE.N5I  >  Y(7,l  ),  YL(  l),nY(l  ) 

IF (T.EO.TCLOiGC  TG  10 
TNTERM  =  EXP(-T) 

TCLD  =  T 
10  CCNTINUE 

R  =  (Y(l,l)  4  Yd, 2)  +  Yd, 3)  +  Yd,4))/2. 

S  =  0. 

DC  20  1=1,4 

20  $  =  S  (R  -  Y(  1,  I)  )*=J‘2/2. 

DC  30  1=1,4 

DY(I)  =  FtNV*Y(2,I)  -  S  +(R  -  Y(1,I))«^2 
CC  25  J=l,4 

25  DY(I)  =  CY(I)  ^  G(I,J)>^Y(1,  J) 

3C  CONTINUE 

CY(5)  =  ♦-INV*(Y(2,5)  +  Y  ( 1 ,  1  )  *Y  (  2 , 6  )  +  Y  (  2  ,  1 )  ^  Y(  1 , 6  )  ) 

CY(6)  =  2.=<'Y(1,6)  +  Y(l,6)^’!«3  -  Yd,l)  *  YL  (1  )  -  l.-T^TER^ 
Cy(7)  =  YL(1)  -  YL(2)  Y(  1,  1)^Y(  1,  6) 

CV(8)  =  YL(1)  -I-  YL(2)  ^  5.^Y(1,1)  *  Y(l,2) 

RETURN 

ENC 


SLBPGUTINE  CERV AL ( Y. YL , T , N, N Y, ^  ,KEPET  ) 
CIMENSICN  Y(7,l  ),  YL(  1)  ,!^(1  ) 

KERcT  =  0 
DG  50  1=1, NY 
5C  Y(2,I)  =  C. 

HINV  =  1. 

CALL  CIFFUN(Y,YL,T,HINV,W) 

CC  100  1=1, NY 
ICC  Y(2,  I)  =  -^d  ) 

RETURN 

END 


36 


Example  2:  This  example  is  a  small  one,  contrived  to  illustrate  the 
possibility  of  derivatives  entering  in  a  nonlinear  fashion.  The  equations 
are 


Yl  -  98y^  +  98y2  =  0 
^2  ■  ®  ^  Yi  +  199y2  =  0 


The  initial  conditions  are 


Yj^  =  Y2  =  1  =  2  . 

Note  that  we  have  supplied  the  explicit  expression  for  the  Jacobian. 
Either  JACMAT  or  a  modified  version  of  DERVAL  must  be  supplied  as  the 
numerical  difference  approximation  to  the  Jacobian  causes  DERVAL  to 
fail. 
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DIMENSICN  Y(7t2) , YL( 1) ,W( 50) 

DATA  N,NY,NL,N,REPS  ,HMAX,HMIi\j,H,T,TENC/3,  2,  1,  2,  1.  E-5 , 25.  ,  1.  E-10, 
1  1 • E “4 f C • # 5 C • / 

Y(l,l)  r  1. 

Y( 1,2)  =  1. 

YL(1)  =  2. 

JSKF  =  0 

CALL  SDESCL ( Y , Y L , T ,T END ,N Y , NL , M, JSKF, 6, 1 , h, HMIN, HM AX, REPS, W ) 
PRINT  6, JSKF 
5  FCRMAT( »0  JSKF  ,14) 

S7CP 

ENC 


SLBROUTINE  CIFFUN  ( Y , YL , f , HINV, 0 Y ) 

DIMENSION  Y(7,1),YL( l),DY(l) 

DATA  T0LC/-79.03/ 

IF(T.EQ.TOLC)GC  TO  10 
TNTERM  =  EXP(“T ) 

TCLD  =  T 
10  CONTINUE 

DY(1)  =  Y(  2,1)’*'HINV  -  9e.*Y(l*l)  +  99.*Y<1,2) 

DY(2)  =  (Y(  2,1)*HINV  )<‘«2  ♦  Y(2,2)*HINV  -  (196.  -  ^  MT  6  RM  )^Y  ( 1 , 1 )  + 
1  199.*Y(1,2) 

DY(3)  ^  VL(  1)  -  Y  (1,1)  -  Y(  1,2) 

RETURN 

ENC 


SL6RDUTINE  JACMAT  (  Y ,  YL  ,  T,  HI  NV,  A2  tN,  NY  ,  E '^5  ,  CY,  FI ,  P  W  ) 
CIMENSTCN  Y(7,1),YL(1),F1(1  ) ,  D  Y  (  1 )  ,  P  ( N  ,  1 ) 

AH  =  A2«HINV 
DG  100  1  =  1, N 
CC  100  J=liN 
100  PVs(  I  ,  J)  ^  6. 

PU(1,1) 

P)^  (1,2) 

PW(2,1) 

PW(2,2) 


-AH 

9S. 


-  98. 


IF(NL.LE.0)RETURN 


-2.*AH^Y( 2,1 )*HINV  -  198.  +  EXP(-r) 
-ah  *  199. 


P^^(3,l) 

PW(3,2) 

PVs(3,3) 

RETURN 

ENC 


1. 
1. 

=  1. 
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Example  3:  This  is  another  contrived  example,  this  one  to  illustrate 
the  use  of  one  type  of  sparse  matrix  storage,  along  with  the  use  of 
iteration  to  solve  the  equations  (2.4).  The  system  of  equations  is 


A  y  +  B  y  =  0  , 


where 


A  = 


f  7  0 

2  8 
0  0 
-3  0 

0  -1 
,  0  0 


0  -3 

0  0 
1  0 
0  5 
0  0 
0  -2 


.3  0  0  .1 

13  0  0 

0  0  10 


B  = 


10 

0 


0 

5 

0 


0  20 
0  0 
0  57 


The  initial  conditions  are 


y^^  =  1000 
y2  =  0 

y3  =  -25 
y4  =  10 
75  =  0 
y^  =  -1000 


0 

0 

0 

0 

4 

0 

0 

0 

0 

0 

6 

0 


0 

0 

0 

0 

6 

-.2 

0 

0 

0 

0 

100 
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The  matrix  storage  scheme  used  for  A,  B,  and  the  Jacobian, 
since  it  has  nonzero  elements  in  the  same  positions  as  A  and  B  ,  is 
that  outlined  by  Gustavson  [4].  Briefly,  one  stores  a  pointer  array 
(here  called  JS)  which  indicates  the  initial  position  of  new  elements 
in  two  other  arrays,  one  of  which  (here  called  JN)  gives  the  column 
number  of  the  element  stored  in  the  corresponding  position  of  the 
coefficient  arrays  (here  called  A  and  B) . 

Thus,  for  the  above  problem  the  arrays  stored  are 


JS: 

1  4 

6  ,^^7^ 

9 _ 

ii__ 

_13 

JN: 

1  4 

6^2 

^  4 

'T~ 

■^6  4 

A  : 

7  -3 

-1  8 

2 

1 

5 

-3 

4  -1 

6  -2 

B  : 

.3  .1 

-.2  3 

1 

1 

20 

10 

6  5 

100  57 

The  elements  of  the 

i— 

row 

are 

Stored  beginning 

at  location 

JS(i) 

of  the 

array  A  , 

and  in  particular 

A(JS(i)  +  k)  is  the  element 

in  the 

1—  row  and  JN(JS(1) 

+  k) 

column  of  A  , 

for  k  =  0 ,  1,  . . 

JS(i+l)  -  JS(i)  -  1  •  For  our  purposes  it  is  necessary  to  access  the 
diagonal  element  easily,  so  we  have  required  that  the  diagonal  element 
be  the  first  element  stored  for  a  given  row.  This  means  that 
JN(JS(i))  =  i  ,  i=l,  ...,  n  .  Note  that  JS(i)  is  the  number  of  nonzero 
elements  in  rows  1  through  i  -  1  , and  that  JS(n+l)  must  be  defined  as 
the  total  number  of  nonzero  elements. 

Problems  similar  to  the  above  arise  when  the  finite  element  method 
is  used  to  discretize  the  space  domain  for  time  dependent  partial  differ¬ 
ential  equations.  Simple  modifications  to  the  subroutines  given  below 
should  permit  solution  of  large  problems  arising  in  that  fashion.  We 
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note,  however,  that  it  is  not  convenient  to  store  synmetric  matrices 
in  this  form  unless  all  nonzero  elements  are  stored.  Storage  of  only 
the  elements  of  the  lower  triangular  matrix  requires  one  to  reference 
columns  of  the  matrix,  which  are  not  readily  accessible.  Even  if  the 
entire  matrix  is  stored,  total  storage  requirements  for  matrices  arising 
in  finite  element  applications  is  still  considerably  less  with  this 
scheme  than  that  required  by  symmetric  band  storage  mode  [5]. 
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ICO 

110 

12C 


3 

4 

5 

6 
7 


C:MMai\‘  /CAT4/-fl(12)  ,B(12),N»  JS(7)  ,JN(12) 

JbtJN 

CATA  T,  TENDfH, JSKF  /  0 .  ,2 30  .  ,  1 . F- 5 , 0  / 

CATA  JSC/1,4,6, 7, c, 11, 13/ 

CATA  JNC/ 1,3, 5, 2, 1,3,4, 1,5, 2, 6, 4/ 

DATA  AP/7#  ,~2*,‘“1»  ,8*,2*,i*,5*,“i«,4»  ,6*,~2*/ 

DATA  BC/.3,.l,-.2,3.,l.,l.,20.,10.,o.,5.,lC0.,57./ 

CATA  YT/1000.,0.,-25. , 10^ , 0. ,“1000./ 

\  =  6 

N  Fl  =  f\  4  1 
JE  =  JSC(NFl)  “  1 
CC  100  1^1, N 
Y(1,I  )  =  Yt  (I  ) 

cc  110  1=1, je 

Ad)  =  AC(I  ) 

ed)  =  ecd) 

JNd  )  =  JNC  d  ) 

C*:  120  1  =  1, NPl 
JS  (I  )  =  JSC  d ) 

PRINT  4,N,  (  JSd  )  .I  =  l,NPl) 

PRINT  5,(Ji\id),I  =  l,  JE) 
point  7, (Ad), 1=1, JE) 

PRINT  6i(e( I ) ,I =1, JE) 

CALL  SDcSOL  (Y,YL,  Td6Nij,N,0  ,N,  JSKP,6,  I ,  H  ,1 .  F-6 ,  1^  5  .  ,  1 .  E  “4  , 

PRINT  3,JSKF 

STCP 

cCRf''AT(  ‘CRETLRN  FRGM  SCcSGL  TH  JSKF  =»,I4) 

FORMAT(»  FOR  THIS  CASE  \=*,I3//»  fHF  JS  ARP  AY  V/ (  I  2 T 10 )  ) 
FGRMAT(/»OTFe  JN  AP R A Y • / / ( I  2  1 10 ) ) 

FCR|^AT(/*0THE  3  ARRAYV/(  12F10.2)  ) 

FORMAT  ( /»0TI-E  A  A  RR  AY  •  /  /  (  1  2  F 10. 2  )  ) 

PNC 


SIBRjUTINC  ClfcjN  ( Y, YL,T, HINV,DY ) 

CCMMON  /CATA/Ad2)  ,8d2),N,  JS(7),  JN(12) 

INTEGER»2  JS,JN 

CIMENSICN  Y  (7, 1 ) , YL(  1) ,0Y( 1 ) 

CC  400  1  =  1, N 
CYd  )  =  C. 

JE  =  JSd) 

JS  =  JS  d  +  1  )  -  1 
DC  300  J=je,JE 

CYd)  =  CY(l)  4  y  (2,  JN(  J)  )=t'A(J)’4'Hli>lV  4  3  (  J  )  «  V  (  1 ,  J  N  (  J  )  ) 
200  CCNTINUE 
400  CCNTINLE 
RETURN 
END 


SLBRGUT  INS  NUITSL  (  PW  ,  DY  ,  F  1  ,  N  , N Y  ,  EP  S  ,  Y Me, X  ,  N W ,  K R  ) 
CCMM3N  /C  AT  A/ Ad2  )  ,8(12  )  ,ND,  JS  (7)  ,  JN(  12  ) 

INTEGER42  JS,JN 

CIMENSICN  Pa( I)  ,CY( I)  , Fl(  I  ), YMAX(  1) 

CATA  CMtG,C^EG^l  /I. 05, .05/ 

KRET  =  0 
SPSS  =  FFS^^2 
eF$A2  =  EPSS’fr.OOOl 
NCIT  =  N 
CC  100  1  =  1, NY 

ICO  Fid)  =  CyI  I  )  /o^i  JS(  I  )  ) 

CC  300  IT=1,NCIT 
RCh  =  0. 

CF  =  0. 

CC  200  1  =  1, NY 
je  =  JS  (I )  4  1 

JE  =  JS(I41)  ^  1 
FN  =  DYd) 

IF(  JB.GT.JE  )GC  Tj  18C 
CC  150  J=JE,JE 

150  FN  =  FN  -  PVs(J)  =«‘F1(  JN(  J)  ) 

130  FN  =  FN'/FW(je-l) 

FN  =  “  F1(I)*GM8GM1 

ACH  =  FI  (I  )  - 


w ) 
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CH  =  CH  4  (  ACri/Y^AX(  I) 

RCH  =  PCH  4  (  ACH/AMAXKABS  (fN),6P^  ) 
200  FUT)  =  FN 

IF (RCH.LT.EPSS)  RETURN 
IF(CH.IT.£PSA2)  RETURN 
300  CCNTINUE 
KRET  =  1 
RETURN 
ENC 


iieROUTINE  JACMAT(Y. YUfT^HINV, A2»N,NY»EPS, CV,F1»PW ) 
CCMMGM  /CATA/A(12)  »6(  12  ) ,  NDUI^,  J  S  (  7  )  ,  JN  ( 12  ) 

INTEGER*2  JStJN 

Clf'ENSICN  Y(7,l)tYL(  DfFK  1),CY<  1),PW(1) 

AH  -  -A2*HINV 
JE  =  JS(N4U  -  1 
CC  100  J=1»JE 

ICO  PU(J)  s  ^H*A(J)  +  B(J) 

RETURN 

ENL 
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Example  4:  This  example  arises  from  a  nonlinear  reactor  dynamics 
problem  where  the  finite  element  method  is  used  to  discretize  the  space 
domain.  The  resulting  system  of  equations  has  the  form 

•  _ 

A  y  -  B  y  +  w(C  y)  y  =  0  , 


where  C  is  a  matrix  with  three  subscripts.  The  i —  equation  can  be 
expressed  as 


N 


N  N 


I 

j=i 


In  this  example 
per  row  in  A 


N  =  28 
and  B  . 


+  w 


I 

3=1 


I 

k=l 


'ijk 


^3 


yv  -  0 


,  and  there  are  at  most  seven  nonzero  elements 
The  nonlinear  term  y .  y,  appears  only  if 

1  K 


both  a^^  and  are  nonzero.  Therefore  a  different  type  of  sparse 

matrix  storage  is  used  for  this  problem. 

An  array,  K  ,  dimensioned  (28,  7)  is  used  to  store  (for  each  row), 
the  columns  subscripts  for  the  nonzero  elements.  For  convenience  in 
accessing  the  diagonal  element,  we  require  that  K(i,l)  =  i  .  We  can 
note  this  matrix  is  simply  the  connectivity  matrix  for  the  finite 
element  grid.  Then  the  nonzero  elements  of  A  and  B  ,  are  stored  in  the 
corresponding  portions  of  the  arrays  A  and  B  respectively.  If  there 
are  in  fact  less  than  seven  nonzero  coefficients  in  a  row,  the  remaining 
K(i,j)  are  set  to  zero. 

The  storage  for  C  is  somewhat  more  complicated.  C  is  symmetric 
(invariant  under  any  permutation  of  subscripts).  The  nonlinear  term  of 


th 

the  i —  equation  was  rewritten  as 
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N  N  _  N  N  _ 

w  I  ^  ‘lik  ^  ^  ^iik  ^k  ’ 

j=l  k=l  j=l  k=j  ^ 


where 


ijk 


‘'ijk 

^Ijk  ^ikj 


k  <  j 
k  =  j 
.  k  >  j 


The  coefficients  then  stored  in  a  (28,28)  array  C  in  the 

order  the  second  and  third  subscripts  are  given  here. 

(K(i,l),  K(i,l)),  ...,  (K(i,l),  K(i,7)),  (K(i,2),  K(i,2)),  (K(i,2), 

K(i,7)) . (K(i,7),  K(i,7))  . 

The  equations  can  then  be  written  in  the  form 


7  7 


-  hj  yKCi.:)’  *  J,  ^KCij)  yu±, 


k) 


=  0 


where 


„  .  k  +  UzlHiizJi  . 

Jk  2 


Because  of  the  large  amount  of  data  for  this  problem  the  input 
arrays  K,  A,  B,  and  C  are  simply  listed  along  with  the  programs  for 
this  example. 
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Clf-ENSICN  Y(7,28)  ,WS(560) 

CCMMON  /CAT4/M28,7)  tB(28,7),C(2e,20)  ,N,NNZ,K(28,7) 

IMEG£«^2  K 

CATA  TPNO,H^IN,HM/5X,EPS,ZOM£GA/  .1,  l.E“12t  .1  ,.01t650.90  3£-U/ 
C/ITA  NN,NY,NL/28, 28tO/ 

^^Z  =  7 

C/iLL  EPRSET  (207,256,-1,1  ) 

C4LL  ePRSET(2C8,256,-l,l) 

N  =  NN 
N  =  N 

NE^C  ==  NNZ+  (NNZ  ♦  l)/2 
PRINT  10 

cc  no  m,NN 

PE4D  1, (K(  1  ,J)  •J  =  1,NNZ) 
no  PRINT  nt(K(I,J)  ,J  =  1,NNZI 
PRINT  12 
CC  120  1=1, NN 
RE/iD  2,  (A(t  ,J)  ,J=1,NNZ) 

V(1,I)  =  0. 

120  PRINT  15,(A(I, J» , JxI.nNZ) 
y(l,l)  =  1.E16 
PRINT  13 
CC  130  !=1,NN 
RE/iC  2,  (e(I,J),J  =  l,NNZ) 

130  PRINT  15,(e(I, J) , J=1,NNZ) 

PRINT  14 
CC  140  1  =  1, NN 
RE^D  2,  (C(  I,J)  ,  J=1,NENDI 
CC  135  J=1,NEN0 
135  C(I,J)  =  d(I,JI  ♦ZOMEGA 
140  PRINT  15,(C(I,J), J=1,NEND) 

JSKF  =  0 
T  =  0. 

H  =  HMIN^IOOO. 

CALL  SDESOL(V,YL,T,TEND,NY,NL,K,  J5KF,6,l,H,h>«IN,H!^AX,  FPS,  ,>£  ) 

PRINT  3, JSKF 

STCP 

1  FCPMAT(16I5) 

2  FORMAT!  7En. 4) 

3  FORMAT ( *0  JSKF  =  • , I3» 

1C  FCRMAT(  ‘IK  ^RPAY« /) 

11  FORMAT (8X, 1418) 

12  FCRMAT(///»0A(I , J )•/) 

13  FORMAT!///* 06(1 ,J)*/) 

14  FC!PMAT!///*CC!I,J)*/) 

15  FCRMAT!8X,1F7E16.6) 

END 


SLBROUTINE  CIFFUN ! Y, YL,T,HINV,DY I 

CCMMON  /CATA/A! 28,7) ,B!28, 7) ,d! 28,28) ,N,NNZ,K!28, 7) 

INTEGERY2  K 

CIMFNSKN  V!7,n,YL!l),DY(l) 

DO  400  1=1, N 
CY!I  )  =0, 

CO  300  J=1,NNZ 

IF (K! I, J).LE.O)GC  TO  310 

QV!I)  =  CY!I)  ♦  Y!2,K!I,J)  )«A(I,J)+HINV  -  B ( I , J ) * Y  !  1  ,  K ( I , J )  )  ♦ 
1  C!I,J)^Y(l,K!l,J))=)'Y!l,K!I,l)i 
30C  CCNTINUt 
310  L  =  NNZ 

CC  360  J1=2,NNZ 
IF!K!I,J1).LE.0)GC  TO  400 
DC  350  J2=J1,NNZ 
L  =  L  +  1 

IF(K! I, J2).LE.0)G0  TO  350 

DY!I)  =  CYd)  +  C!I,L)^Y!1,K(I,J1))>J‘Y!1,K!T,J2)) 

350  CCNTINUB 
360  CCNTINLH 
400  CONTINUE 
RETURN 
END 

SLBPOUTINE  JACMAT  !Y,YL,  T,  HI  NV,A.2,N,  NY  ,  EPS,  CY,F  1,PW  I 
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COMMON  /CAT4/A( 28t  7) tR(28, 7) 28,28)  ,NO,NNZ,K (28, 7) 
INTEGFR^2  K,P 

DIMENSION  Y(7,Uf  YL(1),F1(  I  )  ,  D  Y  (1 )  ,  P  W  (  N  Y  ,  1 ) 
CIMENSUN  P  (7,7) 

DATA  P/AS*0/, lNITP/0/ 

IF  (INITP.EC.NNZ)GG  TQ  99 
IMTP  =  NNZ 
CO  98  L=1,NNZ 
DC  98  M  =  LfNNZ 

98  P(L,M)  =  M  +  (L  -  1)*(2*NNZ  -  L)/2 

99  CCNTINUE 

AH  =-A2«hlNV 
CO  300  1=1, NY 
DC  300  J=1,NNZ 

PWd  ,  J)  =  Ah*A(  !,  J)  -  &(  I  ,  J  ) 

CC  100  1=1, J 

IF(K(I,L).L£.0)G0  TO  295 
ICO  Pi^(I,J)  =  PVn(I,J)  +  C(T,P(L,  J))*Y(1,K(I,L)) 

CO  200  M=J,NNZ 
IF(K( I , M).LE.O)Gu  TO  295 
200  PUdfJ)  =  PW(I,J)  C(I  ,P(  J,M)  )«Y(l,Kd,M)  ) 

295  CCNTINLE 
300  CCNTINUE 
RETURN 
END 


SUBROUTINE  NUI TS  L  ( ,  DY  , '=  1 ,  N ,  NY  ,  E  PS  ,  Y  MA  X  ,  N  E  WP  W ,  KR  E  T  ) 
DIMENSION  PUlNYd  )  ,0Y(1)  ,F1(1  ),YMAX(1  ) 

CCMMCN  /CATfl/A(28,7) ,8(28, 7) ,C (28,28) ,ND,NNZ,K(23, 7) 
CATA  SFC,SPCM1/1 .05,.05/ 

INTEGEP«2  K 
KRET  =  0 
EFSS  =  EFS'<'«2 
EPSA2  =  EPSS'^.OOOl 
NCIT  ^  N 

280  CC  281  1  =  1, NY 

281  Fl( I )  =  CY( I )/PW( 1,1) 

CC  287  IT=1,N01T 

RCF  =  0, 

CF  =  0. 

CO  285  1=1, NY 
FN  =  DY(I) 

CO  28A  J=2,NNZ 

IF(K( I, J).LE.O.OP.K(:, J).GT.NY)GO  TJ  284 
FN  =  FN  -  PVsd,  J)+F1(K(I,  J)  ) 

284  CONTINUE 

FN  =  FN/FU(I,1) 

FN  =  FN*SPb  -  SPDM1*F1( I ) 

ACH  =  FI  (I )  -  FN 

CH  =  CH  4  (  ACH/YMAX(  I  )  )*’!‘2 

PCH  =  RCF  4  (ACH/AMAXl(Ae$(FN),EP$)  )‘^‘’^2 

285  Fid)  =  FN 
IF(RCH.LT.EPSS )Gn  TO  288 
IF(CH.LE.EPSA2)GG  TO  288 

287  CCNTINUE 
KPET  =  3 

288  CCNTINUE 
RETURN 
END 
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INPUT  DATA  FOR  EXAMPLE  4 


16  2 
2  1  6 


4 

3 

e  9 

10  5 

5 

4 

1C  21 

6 

11 

7  2 

1 

7 

2 

6  11 

12  13 

a 

8 

3 

2  7 

13  9 

4 

9 

4 

8  13 

14  15 

10 

10 

5 

4  9 

15  22 

21 

11 

6 

16  12 

7 

12 

7 

11  16 

17  18 

13 

13 

8 

7  12 

18  14 

9 

14 

9 

13  16 

19  20 

15 

15 

10 

9  14 

20  23 

22 

16 

11 

28  17 

12 

17 

12 

16  28 

27  18 

18 

13 

12  17 

27  26 

19  14 

19 

14 

18  26 

25  20 

20 

15 

14  19 

25  24 

23 

21 

5 

1C  22 

22 

21 

10  15 

23 

23 

22 

15  20  ■ 

24 

24 

23 

20  25 

25 

20 

19  26 

24 

26 

19 

16  27 

27 

18 

17  28 

26 

28 

16 

27  17 

8.3776E 

02 

8,3776E  02 

4.18885 

02 

0.0 

0.0 

0.0 

0.0 

5.0265E 

03 

4.1888E  02 

2.09445 

03 

2.51335 

03 

2.09445 

03 

4.18885 

02 

Q.O 

1,6755E 

03 

4.1888E  02 

1.67555 

03 

4.188  85 

02 

0.0 

0.0 

0.0 

5.0265E 

Q3 

4.188eE  02 

2.0944E 

03 

2.51335 

03 

2.09^4E 

03 

4.18885 

02 

0.0 

1.4661E 

03 

4.1888E  02 

1 . 46615 

03 

3.14165 

02 

0.0 

0.0 

0.0 

1.C891E 

04 

2.9322E  03 

4.18885 

03 

2.09445 

03 

8.37765 

02 

0.0 

0.0 

2.84e4E 

04 

2.5133E  03 

4.18885 

03 

6.28325 

03 

6.70215 

03 

6.28325 

03 

4.18d3E 

03 

2. 1782E 

04 

1.6755E  03 

2.0944E 

03 

4.18885 

03 

5. 86435 

03 

4.1888E 

03 

2.0944E 

03 

2.64E4E 

04 

2.5132E  03 

4.18885 

03 

6.28325 

03 

6.70215 

03 

6.28325 

03 

4.18385 

03 

1.9059E 

OA 

1.46615  03 

2.09445 

03 

4.18885 

03 

5. 1313= 

03 

3.1416E 

03 

1.5708E 

03 

2.3457E 

04 

2.9322E  03 

5.02655 

03 

8.37765 

03 

6.28325 

03 

o«o 

0.0 

5.3616E 

04 

6.7021E  03 

8.3776F 

03 

1.04725 

04 

1.08915 

04 

1.04725 

O-^ 

8.37766 

03 

4.6S14E 

04 

5.86435  03 

6.23325 

03 

8.37765 

03 

1.00535 

04 

8.37766 

03 

6.2832E 

03 

5.3616E 

04 

6.70215  03 

8.37765 

03 

1.04725 

04 

1.08915 

04 

1.0^725 

04 

8.3776E 

03 

4.4663E 

04 

5.13135  03 

6.28325 

03 

8.37765 

03 

8.95355 

03 

8.4823E 

0-» 

6.2046E 

03 

3.2515E 

04 

5.02655  03 

5.18365 

03 

1.08125 

04 

1.  04725 

04 

0.0 

0.0 

5.62C8E 

04 

1.08915  04 

1. 08125 

04 

1.19585 

04 

1.19585 

04 

1.08125 

04 

0.0 

8.0582E 

04 

1.00535  04 

1.04725 

04 

1.08125 

04 

1.33125 

04 

1.33125 

04 

2.1284E 

04 

5.62C8E 

04 

1.089  15  04 

1.08125 

04 

1.19585 

04 

1.19585 

04 

1.08125 

04 

0.0 

6.C750E 

04 

8.95355  03 

1.04725 

04 

1.08125 

04 

9.2481E 

03 

1.03485 

04 

9.8764E 

0*3 

3.7699E 

03 

3.14165  02 

1. 57085 

03 

1.88505 

03 

2.9531E 

1.68505  03 

3.14165 

03 

6.2046E 

03 

8.71796 

03 

6.1889E 

04 

8.71795  03 

8.48235 

03 

9. 67645 

03 

1. 24685 

04 

6.3303E 

1.24665  04 

1.03f85 

04 

8.83575 

03 

5.7962E 

04 

9.24815  03 

1.19585 

04 

1.47265 

04 

8986575 

03 

;.4719E 

5^ 

1.19585  04 

1.3312E 

04 

1.76715 

04 

9.4719E 

04 

1.33125  04 

1.19585 

04 

1.47265 

04 

1.76715 

04 

5.3014E 

5.18366  03 

1.47265 

04 

1.1958E 

04 

-1.5816E 

09 

1.17205  09 

1.04495 

09 

0.0 

0.0 

0.0 

0.0 

-3.9622E 

2^ 

1.04495  09 

6.3535E 

08 

4.43385 

09 

3.35355 

08 

1.0^495 

09 

0.0 

-3.1631E 

09 

1.04495  09 

2.34405 

09 

1.04495 

09 

0.0 

0.0 

0.0 

-3.9E22E 

09 

1.04495  09 

6.35355 

08 

4.43385 

09 

6. 35355 

08 

1.04495 

09 

0.0 

-4.3361E 

09 

1.0449E  09 

1.83555 

09 

1.48795 

09 

0.0 

0.0 

0.0 

-6.7925E 

09 

4.56095  09 

6.77785 

09 

6.35355 

08 

1. 17205 

09 

0.0 

0.0 

-1.5222E 

10 

4.4338E  09 

6.77785 

09 

1.90615 

09 

1.12125 

10 

1.90615 

09 

6. 77788 

09 

-1.3585E 

1C 

2.34405  09 

6.35355 

08 

6.77785 

09 

9.12185 

09 

6.77786 

09 

6.3535E 

06 

-i.5223E 

10 

4.4338E  09 

6.77785 

09 

1.90615 

09 

1.12125 

10 

1.90616 

09 

6.7778E 

09 

-2.4104E 

10 

1.83555  09 

6.3535E 

08 

6.77785 

09 

7.33555 

09 

8.44465 

09- 

-6.0318E 

08 

-1.3995E 

10 

4.56095  09 

7.94985 

09 

1.35565 

10 

1.90615 

09 

0.0 

0.0 

-2.9627E 

12 

1.12125  10 

1.35565 

10 

3.17685 

09 

1.79895 

10 

3.17686 

09 

1.5556E 

iO 

-2.7989E 

10 

9.12185  09 

1.90615 

09 

1.35565 

10 

1.59005 

10 

1.35565 

10 

I.9061E 

09 

-2.9627E 

10 

1.12125  10 

1.35565 

10 

3.17685 

09 

1.7989E 

10 

3.17686 

09 

1.3556E 

10 

-4.8828E 

10 

7.33555  09 

1.90615 

09 

1.35565 

1C 

1.02125 

10 

1.1621E 

10 

2.0408E 

09 

-3.5210E 

10 

7.94985  09 

1.36925 

10 

1.60435 

10 

3.17665 

09 

0.0 

0.0 

-8.2408E 

10 

1.7989E  10 

1.60435 

10 

-3.69455 

08 

2.40605 

10 

1.3103E 

1C 

0.0 

-7.2321E 

10 

l.?9Cfl6 

10 

3.  17686 

09 

1.3103E 

10 

1.14766 

10 

1.1476E 

10 

1.62806 

10 

-8.24CaE 

10 

1.7989E 

10 

1.3103E 

10 

2. 4060c 

10- 

•3.69456 

06 

1.6C43L 

10 

0.0 

-e.7ii<;E 

10 

1.0212E 

10 

3. 17686 

09 

1.60436 

10 

2.47986 

10 

3.4658- 

09 

1.33966 

10 

-8.2637E 

09 

1 .48796 

09- 

6.03186 

08 

2.  89536 

09 

-4.7C95E 

10 

2, 89523 

09 

8 .44466 

09 

2.04036 

0  9 

6.4C05E 

08 

-1.0136E 

11 

6.40056 

08 

1.1621F 

10 

1.33986 

10 

4.68226 

09 

-1.3280c 

11 

4.68226 

09 

3.46586 

09 

2.37506 

in 

-  1.  2346E 

11 

2.4796E 

10- 

•3.69456 

08 

7.25336 

09 

2.5750E 

10 

-1. ^2e3E 

1 1 

2.406CE 

10 

1.  14766 

10- 

•3.56886 

09 

-1.4e77E 

11 

1.14766 

10 

2.40606 

10 

3.3929E 

09- 

•3. 56886 

09 

-7. 4644E 

10 

1.36926 

10 

3.39296 

09- 

-3.69456 

08 

4.1888E 

02 

5.58506 

02 

2. 79256 

02 

0.0 

0.0 

0.0 

0.0 

4. 1388E 

02 

2.7925E 

02 

0.0 

0.0 

0.0 

0.0 

1.39636 

05 

0.0 

O.C 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

C.O 

2.3133E 

03 

O.C 

8.37766 

02 

0.0 

5.5850c 

02 

0.0 

0.0 

1.3962E 

02 

2.7925E 

02 

0.0 

0.0 

0.0 

2.7925E 

02 

9.77386 

02 

0.0 

0.0 

0.0 

5.58506 

02 

1.  1 1705= 

5. 5t50E 

02 

0.0 

2.2340E 

C3 

9.7738E 

02 

0.0 

8.37766 

02 

1.39636 

02 

5.58506 

02 

O.C 

9.3776E 

02 

2.79256 

02 

1.11706 

03 

2.7925F 

02 

O.C 

0.0 

J.O 

1.3S63E 

02 

2.7925E 

02 

0.0 

0.0 

0.0 

0.0 

0.37766 

02 

2.7925E 

02 

0.0 

0.0 

0.0 

1.39636 

02 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

2.5123E 

03 

O.C 

8.37766 

02 

0.0 

5.585C6 

02 

C.O 

C.O 

1.3S63E 

02 

2.79253 

02 

0.0 

0.0 

0.0 

2.7925c 

02 

9.7738E 

02 

0.0 

0.0 

0.0 

5.5850E 

05 

1.1170F 

03 

5. 56506 

C2 

C.O 

2.23406 

03 

9.77386 

02 

0.0 

a. 3776E 

02 

1.39636 

02 

5.5650E 

02 

C.O 

4.1888c 

02 

2 .79256 

02 

5.58506 

02 

0.0 

0.0 

0.0 

0.0 

1.3963E 

02 

2.7925E 

02 

0.0 

0.0 

0.0 

0.0 

4.18886 

02 

0.0 

0.0 

C.O 

0.0 

0.0 

0.0 

0.0 

0.0 

O.C 

0.0 

0.0 

0.0 

0.0 

Q.O 

6.70216 

03 

2.2340E 

03 

1.9548E 

03 

8.37766 

02 

0.0 

0.0 

0.0 

1.25666 

03 

1.117CE 

03 

C.O 

0.0 

0.0 

0.0 

1.67556 

03 

5.5850E 

02 

0.0 

0.0 

1.39636 

03 

6. 98136 

02 

0.0 

C.O 

1.1170E 

03 

2.7925E 

02 

0.0 

1.1170F 

03 

0.0 

0.0 

0.0 

1.6755E 

04 

2.2340E 

03 

1.39636 

03 

2.51356 

05 

0.0 

2.2340F 

Od 

1.39626 

02 

8.3776E 

02 

5.58506 

02 

0.0 

0.0 

0.0 

0.0 

1.6755E 

02 

1.11706 

02 

0.0 

0.0 

1.9548E 

03 

2.6529E 

03 

0.0 

0.0 

2.2340E 

02 

2.7925E 

05 

1.39636 

03 

6.42286 

03 

2. 65296 

03 

2. 5133= 

Oi 

5.3058F 

03 

1.3404E 

04 

0.0 

1.11706 

03 

1.95496 

0  3 

4 . 46806 

03 

1.9548E 

03 

8.37766 

02 

5.5650E 

02 

2.79256 

02 

0.0 

0.0 

0.0 

1.9548= 

03 

6.98136 

02 

G.O 

O.C 

0.0 

8.3/766 

02 

1.6755F 

03 

1.1  1706 

0-* 

0.0 

1.95488 

03 

2.51536 

03 

1.11706 

03 

0.0 

1.6755E 

03 

1.95406 

02 

1.81516 

03 

1.6755E 

04 

2.53406 

03 

1.39636 

03 

2.51336 

03 

0.0 

2.2340c 

0^ 

1.3963c 

0  3 

8.37766 

02 

5.585C6 

02 

0.0 

0.0 

0.0 

0.0 

1.67556 

U  3 

1.1170E 

03 

0.0 

0.0 

1.95436 

03 

2.65296 

03 

0  .0 

0.0 

2.23406 

03 

2.7925E 

02 

1.39636 

03 

6.4228E 

03 

2.6529F 

0? 

2. 51356 

03 

5.30536 

03 

6.7C21c 

03 

O.C 

1.11706 

03 

1.95486 

03 

2.23406 

03 

0.0 

0  . .» 

2.79256 

02 

2.7925E 

02 

0.0 

0.0 

0.0 

9.3776E 

02 

6.9813E 

02 

0.0 

O.C 

0.0 

3.37766 

0  2 

1.67556 

03 

1.117CE 

Oj 

0.0 

1.9548E 

03 

1 .5566E 

03 

0.0 

0.0 

0.0 

0.0 

0.0 

1.4242E 

04 

O.C 

3.90956 

03 

3.03036 

03 

2. 5135E 

03 

Q.O 

0.0 

1.  1170E 

03 

O.C 

0.0 

1.11706 

0  3 

O.C 

2.51336 

03 

2.C944F 

03 

1.9548E 

03 

0.0 

0.0 

0.0 

3. 3510F 

03 

1.59636 

03 

0.0 

3.07186 

03 

2.3736E 

03 

0.0 

2.7925F 

03 

O.C 

0.0 

0.0 

3.1835E 

04 

5.565CE 

03 

3.07166 

03 

4.1888E 

03 

0  .  j 

3.90956 

03 

3.07186 

03 

2.5133E 

03 

1.3963E 

03 

0.0 

0.0 

O.C 

0.0 

3.3510E 

03 

1.9548E 

03 

O.C 

0.0 

3.63036 

0  3 

4. 3284E 

03 

0.0 

0.0 

3.90956 

03 

4.46803 

03 

2.23406 

03 

1.06126 

04 

4.3234E 

03 

4.18886 

o:> 

1.03326 

04 

2.84846 

04 

0.0 

2.79256 

03 

3.6303E 

03 

7.8191E 

o3 

3.6303F 

03 

2.5132P 

03 

2.2340E 

03 

1.11706 

03 

0.0 

0.0 

0.0 

6.1436E 

0  J 

2.37366 

02 

0.0 

O.C 

G.O 

2.51336 

0  3 

3.35iOE 

03 

1.95486 

03 

0.0 

4.4680E 

03 

4.19886 

03 

1.95486 

03 

0.0 

3. 25106 

03 

4.46806 

03 

5.1662E 

03 

3.18556 

04 

5.5850E 

03 

3.07186 

03 

4.18886 

03 

0.0 

3.90956 

03 

3.0718E 

C3 

2.5133E 

03 

1.3963E 

03 

0.0 

0.0 

0.0 

O.C 

3.351C6 

03 

1. 95436 

5^ 

O.C 

0.0 

3.63036 

03 

4. 32846 

03 

0.0 

0.0 

3.90956 

03 

4 .468CH 

03 

2.23406 

03 

1.06126 

04 

4.32346 

03 

4.18086 

1.03326 

04 

1.4242E 

04 

O.C 

5.79256 

03 

3.62036 

03 

3.9C95E 

03 

0.0 

0.0 

1.1170E 

03 

1.1170E 

03 

0.0 

0,0 

0.0 

2.51352 

0-> 

2.3736E 

02 

0.0 

O.C 

0.0 

2.5135E 

03 

3.3510'= 

03 

1.95486 

03 

0.0 

4.46806 

05 

5  .09446 

03 

0.0 

0.0 

0.0 

0.0 

0.0 

1.2823E 

04 

O.C 

0.0 

0.0 

4. 18686 

03 

0.0 

O.J 

1.95486 

03 

0.0 

0.0 

1.9548F 

03 

0.0 

4.18886 

08 

0.0 

0.0 

O.C 

0.0 

0.0 

2. 3726F 

03 

2.23406 

03 

C.O 

4.74736 

03 

4.0492E 

03 

0.0 

4.46806 

03 

J.O 

0.0 

0  .0 
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1.A2A2E 

04 

8.93618 

03 

4.7473E 

03 

0.0 

0.0 

4.74736 

03 

0.0 

4* 1868E 

03 

2.2340E 

03 

0.0 

0.0 

0.0 

0.0 

2.3736E 

03 

0.0 

0.0 

g.o 

0.0 

0.0 

O.Q 

0.0 

0.0 

0.0 

0.0 

0.0 

2.3736E 

03 

2.2240E 

03 

0.0 

2.7646E 

04 

0.0 

4.4680E 

03 

0.0 

0.0 

0.0 

4.18806 

03 

3.9095E 

03 

1.9548E 

03 

0.0 

0.0 

0.0 

1.0232E 

04 

4.0492E 

03 

0.0 

O.g 

0.0 

4.18888 

03 

2.3736E 

03 

0.0 

0.0 

6.9813E 

03 

0.0 

0.0 

0.0 

9*9 

0.0 

1.7872E 

04 

1.4242E 

04 

8.9361E 

03 

4.7473E 

03 

0.0 

0.0 

4.7473E 

03 

0.0 

4.  I888E 

02 

2.234CE 

03 

0.0 

0.0 

0.0 

0.0 

2.3736E 

03 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

2.2736E 

03 

2.234CE 

03 

0.0 

1.3823E 

04 

O.C 

4.4680E 

03 

0.0 

0.0 

0.0 

0.0 

1.9548E 

03 

1.9548E 

03 

0.0 

0.0 

0.0 

4.18886 

03 

4.0492E 

03 

0.0 

g.o 

0.0 

4.1888E 

03 

2.3736E 

03 

0.0 

0.0 

6.S812E 

03 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

O.g 

0.0 

0.0 

g.o 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0«0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

g.o 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

g.o 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

3.0 

0.0 

g.o 

O.C 

0.0 

0.0 

0.0 

O.C 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

o.g 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

O.C 

c.o 

0.0 

0.0 

0.0 

0.0 

0.0 

O.C 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

O.C 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

O.C 

0.0 

0.0 

o.g 

0.0 

0.0 

0.0 

0.0 

0.0 

3.0 

O.C 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

3.0 

0.0 

3.0 

0.0 

50 


References 


1.  R.  L.  Brown  and  C.  W.  Gear,  "Documentation  for  DFASUB  -  a  program 
for  the  solution  of  simultaneous  implicit  differential  and  nonlinear 
equations,"  Report  no.  UIUCDCS-R-73-575 ,  University  of  Illinois  at 
Urbana  -  Champaign,  Urbana,  Illinois,  July  1973. 

2.  C.  W.  Gear  Numerical  Initial  Value  Problems  in  Ordinary  Differential 
Equations ,  Prentice-Hall,  Englewood  Cliffs,  NJ,  1971. 

3.  C.  W.  Gear,  "Simultaneous  Numerical  Solution  of  Differential-Algebraic 
Systems,"  IEEE  Trans,  on  Circuit  Theory,  CT-18(197l)89-95. 

4.  F.  G.  Gustavson,  "Some  Basic  Techniques  for  Solving  Sparse  Systems 

of  Linear  Equations,"  pp  41-52  in  Sparse  Matrices  and  Their  Applications, 
D.  J.  Rose  and  R.  A.  Willoughby,  eds..  Plenum  Press,  New  York  - 
London,  1972. 

5.  D.  Salinas,  D.  H.  Nguyen,  R.  Olsen,  and  R.  Franke  "An  Optimal  Compact 
Storage  Scheme  for  Nonlinear  Reactor  Problems  by  FEM,"  Manuscript 

(to  be  submitted) . 


51 


Distribution  List 


No,  of  copies 

Defense  Documentation  Center  12 

Cameron  Station 
Alexandria,  VA  22314 

Library  2 

Naval  Postgraduate  School 
Monterey,  CA  93940 

Dean  of  Research  2 

Naval  Postgraduate  School 
Monterey,  CA  93940 

Naval  Postgraduate  School 
Department  of  Mathematics 

Ladis  D.  Kovach,  Chairman  1 

Professor  C.  Comstock  1 

Professor  F.  Faulkner  1 

Professor  R.  Franke  10 

Naval  Postgraduate  School 
Department  of  Mechanical  Engineering 

Professor  D.  Salinas  1 

Professor  D.  Nguyen  1 

Professor  R.  Newton  1 

Professor  G.  Cantin  1 

Naval  Postgraduate  School 
Department  of  Aeronautics 

Professor  D.  Collins  1 

Professor  R.  Ball  1 

Naval  Postgraduate  School 
Computer  Center 
Monterey,  CA  93940 

Professor  D.  Williams  1 

Mr.  Roger  Hilleary  1 

Dr.  Richard  Lau  1 

Office  of  Naval  Research 
Pasadena,  CA  91100 

Chief  of  Naval  Research  2 

ATTN:  Mathematics  Program 
Arlington,  VA  22217 


52 


Mr.  W.  J.  Dejka,  Code  4000 

Naval  Electronics  Laboratory  Center 

San  Diego,  CA  92152 


1 


Argonne  National  Laboratory 

Argonne,  IL  60439 

ATTN:  Mr.  Gary  Leaf,  AMD 


1 

1 


Mr.  Tilak  Chawla,  RAS 


Professor  C.  W.  Gear 
University  of  Illinois 
Urbana,  IL  61801 


1 


Dr.  A.  G.  Hindmarsh 
Lawrence  Livermore  Laboratory 
Livermore,  GA  94550 


1 


Sandia  Laboratories 


Albuquerque,  NM  87115 
ATTN:  D.  A.  Dahlgren 


1 

1 


L.  F.  Shampine 


Mr.  R.  E.  Huddleston 


1 


Sandia  Laboratories 
Livermore,  CA  94550 

Air  Force  Weapons  Laboratory 
Kirtland  AFB 
Albuquerque,  NM  87115 

ATTN;  C.  M.  Walters  1 

CAPT.  C.  W.  Stein 

Mr.  R.  D.  Birkhoff  1 

Oak  Ridge  National  Laboratory 

Union  Carbide  Corporation 

P.  0.  Box  X 

Oak  Ridge,  TN  37830 

Mr.  R.  M.  Sternheimer  1 

Brookhaven  National  Laboratory 
Upton,  NY  11973 

Mr.  B.  F.  Maskewitz  1 

RSIC,  Neutron  Physics  Division 
Oak  Ridge  National  Laboratory 
Oak  Ridge,  TN  37831 

Mr.  J.  L.  Black  1 

U.  S.  Naval  Research  Laboratory 
Code  770 

Washington,  DC  20390 


53 


Los  Alamos  Scientific  Laboratory 

P.  0.  Box  1663 

Los  Alamos,  NM  87544 

ATTN;  S.  Evans,  T-6  1 

C.  Young,  J-14  1 

Professor  R.  E.  Barnhill  1 

Department  of  Mathematics 

University  of  Utah 

Salt  Lake  City,  Utah  84112 


LT  Donald  Kinsman  1 

Fleet  Numerical  Weather  Central 
Monterey,  CA  93940 


54 


U 17 547 5 


DUDLEY  KNOX  LIBRARY  -  RESEARCH  REPORTS 


5  6853  01071186  4 


